Theoretical insights into electronic nematic order, bond-charge orders, and plasmons in cuprate superconductors

In this article, we focus on the charge degree of freedom in cuprate superconductors and review theoretical insights into the electronic nematic order, bond-charge orders, and plasmons. The low-energy charge dynamics is controlled by the spin-spin interaction J, which generates various bond-charge ordering tendencies including the electronic nematic order. The nematic order is driven by a d-wave Pomeranchuk instability and is pronounced in the underdoped region as well as around van Hove filling in the hole-doped case; the nematic tendency is weak in the electron-doped region. Nematicity consistent with the d-wave Pomeranchuk instability was reported for hole-doped cuprates in various experiments. Although the t-J and Hubbard models correctly predicted the proximity to the nematic instability in cuprates far before the experimental indications were obtained, full understanding of the charge ordering tendencies in hole-doped cuprates still requires further theoretical studies. In electron-doped cuprates, on the other hand, the d-wave bond-charge excitations around momentum q=(0.5pi,0) explain the resonant x-ray scattering data very well. Plasmon excitations are also present and the agreement between the large-N theory of the t-J model and resonant inelastic x-ray scattering measurements is nearly quantitative in both hole- and electron-doped cuprates. Theoretically the charge dynamics in cuprates is summarized as a dual structure in energy space: the low-energy region scaled by J, where the nematic and various bond-charge orders are relevant, and the high-energy region typically larger than J, where plasmons are predominant.

where N s is the total number of lattice sites, S i = 1 2c † iα σ αβciβ with Pauli matrices σ, andc † iσ (c iσ ) is the creation (annihilation) operator of electrons with spin σ at site i in the restricted Hilbert space where the double occupancy of electrons is prohibited at any site. After the Fourier transformation, we obtain Here = 1 2 (cos k x + cos k y )(cos k ′ x + cos k ′ y ) implying four different channels.
Let us focus on the forward scattering processes and put q = 0 in Eq. (2); see Sec. 3 for a general q 0. In this case, the first term in Eq. (4) describes the s-wave channel and corresponds to the so-called uniform resonating-valence-bond (RVB) order in the t-J model. 45,46) The second term is the d-wave channel and describes the nematic order. The importance of this d-wave channel was not recognized until 2000. 42,43) The third and fourth terms are not relevant in the presence of inversion symmetry. Therefore the nematic interaction can be extracted from the spin-spin interaction as where d k = cos k x − cos k y . The nematic order parameter may be defined as When the antiferromagnetic interaction (J > 0) is considered, the nematic channel becomes attractive. The functional form of the right-hand side of Eq. (5) is the same as the so-called Landau interaction in the d-wave spin-symmetric channel if we identifyc kσ with the usual electron operator c kσ . In this sense, the electronic nematic instability described by Eq. (5) is often called as a d-wave Pomeranchuk instability (dPI), referring to his paper 47) about the stability (not instability) condition of Fermi liquids. Note that the nematic order can occur even without breaking his stability condition because the transition can be of first order at low temperature. 48,49) Equation (5) is not a special feature of the spin-spin interaction. It is easy to see that the nearest-neighbor Coulomb interaction (V > 0) on a square lattice contains the same If one starts with the Hubbard model, the interaction on the right-hand side of Eqs. (5) and (7) is generated by electron-electron interactions as a low-energy effective one. 44) One can also obtain the dPI in a continuum model with central forces. 51,52) Given that the functional form of the Landau interaction is quite general and may describe the forward scattering interaction in different electron systems, the presence of the d-wave Pomeranchuk interaction itself is quite general, independent of whether the system is defined in the strong coupling limit as in the t-J model or in a weak coupling model.

Typical phase diagram
Obviously the spin-spin and Coulomb interactions contain other ordering tendencies, too and the nematic order is regarded as one of them. But first let us focus on the nematic order and clarify its typical property.  Figure 1 is the phase diagram obtained in the slave-boson mean-field theory in the t-J model by discarding orders competing with the nematic order. The dPI occurs in the lowdoping region and is pronounced upon approaching half-filling. The energy gain is described by 43) and Here f (x) is the Fermi distribution function and ξ k is the renormalized electron dispersion.
That is, the nematic order occurs when the a term exceeds unity. This a term describes the d-wave weighted density of states averaged over an energy interval of order of temperature T around the chemical potential and becomes large in two different cases. One case is that the band width becomes narrower. In the slave-boson mean-field theory, the nearest-neighbor hopping integral t is renormalized to be tδ, where δ is doping rate measured from half-filling.
The renormalization of t to tδ is a special feature of the strong electron correlations that the double occupancy of electrons is prohibited at any lattice site. This is the major reason why the onset temperature increases with decreasing doping in Fig. 1. The other case is that the system is close to van Hove filling. In Fig. 1 the van Hove filling is located around δ = 0.10.
The enhancement of the nematic instability there is not visible, implying that the effect of the band narrowing is dominant. Similar results to Fig. 1 were also obtained in the strong coupling Hubbard model. [53][54][55] For different choices of band parameters, the phase diagram does not change qualitatively close to half-filling. 43) An additional feature is that van Hove filling can be located in a large doping region and the nematic order occurs also around the van Hove filling, but with temperature much smaller than Fig. 1 (Refs. 37,43). The phase digram becomes similar to that obtained in a weak coupling model except for the temperature scale (see Fig. 2

below).
In a weak coupling model, the nematic order occurs around van Hove filling with a domeshaped transition line as shown in Fig. 2: a second-order transition on the roof and a first-order transition near the edges of the dome. 48,49) The end points of the second-order transition are tricritical points.
In the weak coupling limit, the phase diagram Fig. 2 is fully determined by a single energy scale. 49) For example, the transition temperature at van Hove filling is given by where γ = 0.577 is the Euler constant,ḡ is the dimensionless coupling constant, and ǫ Λ corresponds to the typical energy scale around the saddle point contributing to the nematic   (7)] with the kinetic energy ξ 0 k = −2t(cos k x + cos k y ) − 4t ′ cos k x cos ky − µ with t ′ /t = −1/6 and the chemical potential µ. T 2nd c is a second-order transition line at high temperature and T 1st c describes first-order transition lines at low temperature. The end points of the second-order transition line are tricritical points. "T 2nd c " is fictitious second-order transition lines preempted by the first-order transition. µ = −2/3t (dotted line) corresponds to van Hove filling. Adapted from Ref. 49. instability. The tricritical temperature is with α = 0.4515. Hence the dimensionless ratios of different quantities become universal: Note that the functional form of Eq. (10) is exactly the same as T c in the BCS theory 56) and the presence of universal ratios are also shared with the BCS theory.
While a certain approximation and some simplification are usually needed to compute the phase diagram, we emphasize that the presence of the dPI interaction in the attractive channel itself does not depend on an approximation. Higher order corrections to the approximation may modify quantitative features of the phase diagram, but may not introduce a drastic change as long as the nematic instability survives. In fact, exact diagonalization 57) and variational Monte Carlo 58) studies of the t-J model suggest qualitatively the same phase diagram as of the presence of the attractive interaction. 60) Note that the nematic instability is associated with the breaking of Ising symmetry and thus it can occur even at finite temperatures in the two dimensions.

Competition with other ordering tendencies
The electron-electron interaction can also drive other orders such as superconductivity, antiferromagnetism, charge-density-wave, and bond-charge orders. Which order can be the leading one in the t-J and Hubbard models? In the t-J model, d-wave pairing is a stronger instability than the nematic order in the slave-boson mean-field theory 43) and the variational Mote Carlo study. 58) Antiferromagnetism also preempts the nematic instability in the slaveboson mean-field theory. Exact diagonalization 57) found a strong tendency of the nematic instability, but could not conclude that the nematic order is indeed the leading instability.
The nematic instability is usually weaker than bond-charge orders in a large-N theory, 37,61) but can become a leading one in a heavily doped region when a large |t ′ | is introduced. 37) In the strong coupling Hubbard model, a strong nematic tendency was found close to a Mott transition in the cellular dynamical mean-field theory 53) and the dynamical cluster approximations. 53,54) Spontaneous symmetry breaking to the nematic phase was obtained by considering a coupling to the lattice in the cellular dynamical mean-field theory. 55) In the weak coupling Hubbard model, extensive calculations in different functional renormalization group schemes showed that the d-wave superconductivity and antiferromagnetism are the leading ones [62][63][64][65] and that charge-density-wave becomes a leading order when the sizable nearest-neighbor Coulomb repulsion is taken into account. 65) On the other hand, the coexistence of the nematic order and d-wave superconductivity was found in the second-order perturbation theory 66) and the dynamical mean-field theory combined with the fluctuation exchange method 67) in the Hubbard model. Competition of the nematic instability and superconductivity was studied 68) in a model, which contains only the BCS pairing interaction and the forward scattering interaction such as Eqs. (5) and (7).

Big response to a small anisotropy
Available theoretical results imply that the nematic tendency is surely present in both t-J and Hubbard models. The point here is that even if the nematic tendency is preempted by a different order, the susceptibility of the nematic order can remain large. In this case, the electronic property becomes very susceptible to a small external anisotropy and exhibits a big anisotropy.
The big response to a small anisotropy is a robust and a general property of the nematic physics when the system is close to the nematic instability. This physics was extensively studied in the slave-boson mean-field theory of the t-J model 42,43,[69][70][71][72] and in the Hubbard model in different schemes: the cellular dynamical mean-field theory 53,55) and the dynamical cluster approximation. 53,54) Among various results, we here present Fig. 3(a), which shows the degree of the anisotropy as a function of temperature T for various choices of doping δ in the t-J model with 3 % anisotropy in the nearest-neighbor hopping integral between the x and y directions. 43) At there is a small anisotropy coming from the original external anisotropy. With decreasing temperature, the anisotropy is strongly enhanced because of the underlying nematic correlations and takes a cusp at the onset temperature of the d-wave pairing gap. The competition with the pairing formation then suppresses the anisotropy. Nevertheless, the big anisotropy remains at zero temperature, which is a few times larger than that at high temperature. This enhancement is pronounced for lower doping, which is easily expected from  Fig. 3(b). Since the external anisotropy is rather small in Fig. 3(a), the curves in Fig. 3(a) are regarded as the temperature dependence of the nematic susceptibility.
In the theoretical scheme in Ref. 43, the nematic tendency increases with decreasing doping as shown in Figs. 1 and 3(a). However, results in the vicinity of δ = 0 should be taken carefully because the Mott physics is considered mainly as the band narrowing effect in the slave-boson mean-field theory of the t-J model. Emergence of antiferromagnetism is not considered near half-filling. As described in Sec. 2.3, the nematic order competes with other ordering tendencies and thus would be suppressed eventually in the vicinity of half-filling.
The strong coupling Hubbard model with a large on-site Coulomb interaction U (U/t = However, it should be kept in mind that the possible magnetic instability was not considered in those studies.

Spectral function
Hereafter we do not distinguish between electron operatorsc kσ and c kσ because the nematic physics is relevant to both strong and weak coupling models. The nematic correlation function is defined as where q and ω are momentum and energy transfer, respectively, Γ is a positive infinites-  reported in cuprates despite many experimental indications of nematic correlations. 21,[74][75][76][77][78][79][80][81][82] The linear dispersion shown in the inset in Fig. 4(a) should not be confused with the zero-sound mode, which is also characterized by a gapless linear dispersion near q = (0, 0).
The zero-sound mode is driven by a short-range repulsive interaction and is realized above the particle-hole continuum. 83) On the other hand, the nematic mode originates from a shortrange attractive interaction in a d-wave channel [see Eqs. (2) and (4)] and is realized inside the particle-hole continuum in the normal phase.

Electron self-energy from nematic fluctuations
In the normal phase, the nematic fluctuations appear in a low-energy region as shown in Fig. 4(a). Its propagator can take the form in general 84,85) where ν n = 2πnT is a bosonic Matsubara frequency; ξ is the nematic correlation length, ξ 0 and u are non-universal parameters determined by the momentum dependence of the interaction strength and the band structure; g is a coupling strength at q = (0, 0). Quantum nematic fluctuations are described by finite Matsubara frequencies ν n 0 in Eq. (14). References 84 and 86 studied how they renormalized the one-particle property of electrons. Non-Fermi liquid behavior was obtained at a quantum critical point, otherwise the Fermi liquid state was stabilized in the ground state. In both cases, the spectral function for one-particle excitations A(k, ω) exhibits a single peak at k = k F and ω = 0; k F is the Fermi momentum.
At finite temperature, a more drastic effect occurs. 85) Since the quantum fluctuations are cut off by temperature, we may focus on the thermal fluctuations in a relatively high temperature region and put ν n = 0 in Eq. (14).
The self-energy computed perturbatively to first order is shown in Fig. 5(a) for several choices of the nematic correlation length ξ. The self-energy exhibits a peak at ω = 0 for the momentum on the Fermi surface and the peak develops to be sharper with increasing ξ. Consequently, as shown in Fig. 5(b), a peak of the spectral function A(k F , ω) at ω = 0 splits and forms a double peak structure with strong suppression of the spectral weight at ω = 0, indicating a gaplike feature around the Fermi surface, although the system is in the disordered phase. We refer to this gaplike feature as a pseudogap, which is similar to the pseudogap observed in cuprates. 24) The results in Fig. 5 suggest a pseudogap driven by the thermal nematic fluctuations.
However, when more precise calculations were performed by including the self-consistency in the perturbative analysis and also vertex corrections, 85) the pseudogap feature in Fig. 5(b) disappears. The obtained self-energy is shown in Fig. 6(a). It still exhibits a peak at ω = 0 for k = k F and the peak is pronounced more upon increasing ξ. But compared with Fig. 5(a), the absolute value of ImΣ(k F , ω) is suppressed substantially. Consequently the spectral function forms a single peak at ω = 0 as shown in Fig. 6(b), where there is no indication of the gaplike feature with increasing ξ, although the peak at ω = 0 is suppressed.
In Fig. 6(b) the peak width is proportional to log ξ and is broadened with increasing ξ. This feature is very different from that from the quantum nematic fluctuations. In the quantum critical regime, 84) the spectral function also exhibits a single peak, but with the peak width proportional to T ξ; here ξ diverges as (T | log T |) −1/2 upon approaching the quantum critical point at T = 0. The single peak structure in Fig. 6(b) should not be understood as the indication of a Fermi liquid. Rather it can be interpreted as an incoherent peak in the sense that the self-energy exhibits a peak at ω = 0 as shown in Fig. 6(a).   While the nematic correlation length ξ is rather small, the spectral weight is substantially suppressed near k = (π, 0) and exhibits a Fermi-arc-like feature, although no gaplike feature is realized around k = (π, 0) as seen in Fig. 6(b).

Probes of nematicity 2.7.1 Raman scattering and ultrasound
Electronic Raman scattering in the B 1g symmetry measures directly the nematic correlation function for q = 0 in Eq. (13) and can provide microscopic evidence of nematic fluctuations in real materials. 87) When the system approaches the nematic instability with decreasing temperature, the so-called central peak was predicted [ Fig. 7(a)]. On the other hand, when the nematic quantum critical point is hidden inside the superconducting phase and the system approaches it upon changing the doping rate, the softening of the B 1g peak was predicted [ Fig. 7(b)].    (7)] and also with phonons; the electron self-energy was also considered in the normal phase by modeling the bosonic spectral function α 2 F phenomenologically. B 1g phonons also couple directly to the nematic fluctuations and Raman scattering exhibits a characteristic feature near the nematic instability. 87) A caveat here is the B 1g phonon energy is assumed to be relatively large, e.g., around 40 meV in YBa 2 Cu 3 O 6+y (YBCO). In this case, the original phonon frequency itself does not change much and stays around ω ≈ 0.25 in Fig. 7(c) and 0.28 in Fig. 7(d). Instead, a central peak emerges in the phonon spectrum in the normal phase [ Fig. 7(c)] and a soft phonon mode in the superconducting phase [ Fig. 7(d)].
When the B 1g phonon's frequency is very small, the predicted double peak structure in the normal phase in Fig. 7(c) may overlap with each other and look like a single peak. In the superconducting state, the original phonon mode simply softens down to zero energy and no additional low-energy peak emerges.
In Bi 2 Sr 2 CaCu 2 O 8+δ (Bi2212) around van Hove filling, where the pseudogap temperature nearly vanishes, the B 1g Raman scattering 81) showed a central peak as predicted in Fig. 7(a).
The nematic correlations are then suppressed inside the pseudogap region. Further exploration by Raman scattering for different cuprate compounds will serve to elucidate how the nematic fluctuations evolve in the cuprate phase diagram. 6) The theory of Raman scattering from nematic fluctuations can be applied to other materials. 88 93) and frustration of localized spins. 94) Ultrasound waves also couple to the nematic fluctuations. In particular, it was pointed out theoretically that the nematic fluctuations enhance the transverse sound attenuation and sound-velocity softening along the [110] direction upon approaching the nematic instability whereas they remain unaffected along the [100] direction. 95) The transverse phonons along the [110] direction have the electron-phonon vertex characterized by B 1g symmetry, the same symmetry as the B 1g phonon discussed above. We are not aware of experimental tests of the ultrasound wave anomaly in cuprates.

xy-anisotropy of physical quantities
YBCO has an intrinsic small anisotropy coming from the orthorhombic crystal structure.
Various physical quantities, however, can exhibit a large xy-anisotropy due to a coupling to the underlying nematicity, which is a robust feature of the nematic physics as was discussed in Fig. 3.
The first experimental signature was obtained in the xy-anisotropy of the magnetic excitation spectra in YBCO. 21,[74][75][76] While the anisotropy is observed in the magnetic excitations, this does not necessarily mean that the origin of the nematicity should lie in the magnetic sectors such as spin nematic. [96][97][98] In fact, the observed anisotropy was well captured in terms of the underlying nematic correlations from the dPI in the t-J model. 70,72) The Nernst coefficient also exhibits the strong xy-anisotropy, which seems pronounced below the pseudogap temperature 77) or below a certain temperature inside the pseudogap phase. 78) The large anisotropy of the Nernst coefficient can be undertood as a consequence of Fermi surface distortions due to nematicity. 99) Magnetic torque measurements for YBCO reported that a component of the two-fold symmetry starts to enhance across the pseudogap temperature and the authors argued that the pseudogap temperature corresponds to the nematic instability. 79) If this is indeed so, the electronic Raman scattering in B 1g symmetry should exhibit a central peak close to the pseudogap temperature [see Fig. 7(a)]. On the other hand, in contrast to Fig. 5(b), the elaborate calculations in Fig. 6 found that the nematic fluctuations do not generate a gap along the FS, although the Fermi-arc-like feature is generated. In addition, the experimental data 79) can also be interpreted in terms of antiferromagnetic fluctuations with the Dzyaloshinskii-Moriya interaction without invoking the nematic physics. 100) The magnetic torque results will stimulate further studies on the relationship between the nematic order and the pseudogap.

Global view of the nematic physics: a concept of Griffiths wings.
The nematic order couples directly to an external xy anisotropy such as strain, uniaxial pressure, and an orthorhombic crystal structure. In the presence of xy anisotropy, the nematic order parameter becomes finite and thus no second-order phase transition occurs. However, a first-order phase transition is still possible. The nematic phase diagram was revealed in the three-dimensional space spanned by the chemical potential µ, the external xy anisotropy µ d , and temperature T (Ref. 101). The inset in Fig. 8 is a mean-field phase diagram. The phase diagram is symmetric with respect to the axis µ d = 0 and almost symmetric with respect to the axis of µ = 0 as long as |t ′ | is small. A wing (colored in orange) develops from the first-order phase transition line at µ d = 0. Crossing this wing, the nematicity changes discontinuously, that is, a meta-nematic transition occurs. This is a first-order nematic phase transition. The  When weak nematic order-parameter fluctuations are included, the wing splits into two as shown in the main panel in Fig. 8. One wing is realized in a large anisotropy region and the other is in a region close to zero anisotropy. In this case, a QCEP is realized close to the tetragonal phase (µ d = 0) and thus strong nematic fluctuations are expected even if the external anisotropy is present. No phase transition is present between the two wings depicted by the dashed line in Fig. 8.
At zero temperature, YBCO is expected to locate near the point of QCEP 2 in Fig. 8 and µ = 0 may correspond to half-filling. In fact, the magnetic excitation spectra exhibit the pronounced anisotropy in YBCO 6.3 , YBCO 6 Fig. 8 corresponds to a carrier density between YBCO 6.45 and YBCO 6.6 .
Referring to a pioneering work by Griffiths about the wing structure associated with a first-order phase transition, namely a concept of tricritical point, in the He 3 -He 4 mixtures, 102) the wing structure in Fig. 8

Bond-charge orders and their fluctuations
The possible importance of bond-charge orders in cuprates was already recognized in early theoretical studies. [105][106][107][108][109] The RVB theory of the t-J model was formulated also by introducing bond-order parameters. 45,46) Recent resonant x-ray scattering (RXS) 18,110) and resonant inelastic x-ray scattering (RIXS) 111) measurements suggest a bond-charge order, which however seems different from what was discussed in early studies.
Extensive studies of possible charge orders were performed in a large-N theory of the t-J model on a square lattice at leading order. 37,61) This theory has an advantage to study all possible charge instabilities present in the t-J model on an equal footing. The theory can be formulated in different schemes 112,113) including a path integral representation. 114) The path integral formalism was shown to yield results consistent with exact diagonalization. 115) We review theoretical insights obtained in that formalism focusing on bond-charge orders in this section. Usual on-site charge fluctuations shall be reviewed in the next section by including the long-range Coulomb interaction.
In the leading-order theory of the large-N expansion, the number of spins is extended from D ab (q, ω) are obtained as a 6 × 6 matrix and given by where a and b run from 1 to 6; q and ω are momentum and energy transfer, respectively. The quantity D (0) ab describes bare charge susceptibilities and is renormalized by the boson selfenergies Π ab at the order of 1/N. Explicit forms of D (0) ab and Π ab are given in Ref. 37. In this scheme, the tendency toward phase separation becomes rather strong and thus the nearestneighbor Coulomb interaction is usually introduced to avoid the phase separation. This does not introduce essential changes in the underlying tendency of various bond orders as well as their excitation spectra.

Bond-charge orders
The instability of the paramagnetic phase is signaled by the divergence of the static susceptibilities D ab (q, 0) for a continuous phase transition. Therefore we study eigenvalues and eigenvectors of the matrix [D ab (q, 0)] −1 . When an eigenvalue crosses zero at a given doping rate, temperature, and momentum, the instability occurs toward a phase characterized by the corresponding eigenvector. Among numerous possibilities, essentially only a few bondcharge instabilities are relevant to parameters appropriate to cuprates. i) Flux phase with q ≈ (π, π). In this phase, currents flow in each plaquette; see Fig. 9(a). ii) the dPI [see Fig. 9(b)]. The dPI leads to the electronic nematic state as already described in Sec. 2. The dPI can have a finite momentum q near (0, 0), which is often referred to as an incommensurate dPI. 65,[116][117][118][119] iii) Various bond-charge orders with q ≈ (π, π): uniaxial bond-charge order along the x or y direction, whereas sbond and dbond with q = (π, π) have a bond amplitude modulated along both x and y directions, and its amplitude is inphase and antiphase, respectively. The dbond with q = (0, 0) is the same as the dPI. That is, the dPI is a special type of the dbond among various types of bond orders.
The phase diagram obtained in the large-N theory of the t-J model 61) is shown in Fig. 10 for both hole-and electron-doped cases. At half-filling, flux phase, dPI, dbond, sbond, and unibond have the same onset temperature T c = J/8. Upon carrier doping such a degeneracy is lifted. On the hole-doped side, flux phase with q = (π, π) is the leading instability and the dPI with q ≈ (0, 0) is the second leading one; when a larger |t ′ | is taken, the dPI would extend to a higher doping region and become the leading one there. 37) The sbond and unibond have ordering tendencies around q = (π, π), which are nearly degenerate and much weaker than flux phase and the dPI.
The charge ordering tendency exhibits a strong particle-hole asymmetry. In contrast to the hole-doped case, the charge ordering tendency is limited only close to half-filling in the electron doping region. The dPI becomes the weakest instability and flux phase is leading at high temperature. At low temperature in a moderate doping region, dbond with q ≈ (0.8π, 0.8π) becomes dominant and the ordering tendency of sbond and unibond with q ≈ (π, π) are located close to the dbond.
In the large-N theory, the effect of spin fluctuations does not appear at the leading order and thus no magnetic instability is observed in Fig. 10. While this is an advantage of the large-N theory in that one can focus on the charge degree of freedom in the t-J model, a comparison with experiments should be made carefully by keeping in mind a possible magnetic instability.
On the hole-doped side, one may assume the magnetic order close to half-filling, for example, in δ 0.05. Hence the theory predicts the flux instability in a large doping region, which is however not observed in experiemnts. This inconsistency remains to be studied.
Since the second leading one is the dPI, we can expect that the systems has a large nematic susceptibility and becomes sensitive to an external xy anisotropy even if the dPI does not occurs. This feature is in fact consistent with the nematic tendency observed in cuprates; 21,[74][75][76][77][78][79][80][81][82] see also Fig. 3.
On the electron-doped side, the magnetic phase may extend to a wide doping region, which could hide all charge ordering tendencies or make only dbond, sbond, and unibond relevant to the reality. As we shall show below, the dbond ordering tendency is consistent with the RXS 18,110) and RIXS 111) measurements.

Bond-charge excitations
The charge excitation spectrum is computed from Eq. (15). To extract the bond-charge component, we need to project D −1 ab onto the corresponding eigenvector. We obtain the dbond, sbond and flux order susceptibilities as follows: The difference between χ dbond and χ sbond lies in the sign in front of D −1 34 and both quantities become identical for q = (π, q y ) and (q x , π). Since Fig. 10 cannot capture the charge order around q = (0.6π, 0) reported in RIXS for the hole-doped cuprates, [15][16][17] we here present theoretical results for the electron-doped cuprates. Theoretically the electron-doped cuprates are addressed by taking a positive value of t ′ /t in the t-J model. 120,121) Bond-charge excitation spectra of χ dbond (q, ω), χ sbond (q, ω), and χ flux (q, ω) are shown in Fig. 11 along the symmetry axes. χ dbond exhibits large spectral weight at low energy around q = 0.8(π, π). This spectral weight is associated with the leading dbond instability at δ c = 0.133 in Fig. 10. Along the direction (0, 0)-(π, 0) the spectrum has rather high intensity and its energy goes down toward the momentum q = (0.5π, 0). χ sbond shows a low-energy dispersion around q = (π, π), which is related to the proximity to the corresponding instability at δ c = 0.114 (see Fig. 10). Its spectral weight disperses upwards forming a V-shape and loses intensity with increasing ω. In contrast to the case of χ dbond , there is no ordering tendency along the (0, 0)-(π, 0) direction.
The important insight obtained here is that the low-energy charge excitations are not dominated by a certain bond-charge order, but by various types of bond-charge orders especially around q = (π, π). RIXS measurements in such a region have not been performed. What was extensively studied is a region along the (0, 0)-(π, 0) direction by RXS. 18,110) While a signature of flux phase has not been reported, RXS 18,110) and RIXS 111) experiments reported the charge ordering tendency around q = (0.5π, 0) as seen in Fig. 11(a). Detailed theoretical studies 123,124) showed that the temperature and doping dependences of the spectrum are well captured in terms of the dbond charge excitations.
By taking the layered structure in cuprates and the long-range Coulomb interaction into account, the large-N theory of the layered t-J model on a square lattice turned out to explain the charge excitations observed around q = (0, 0) in terms of acoustic-like plasmons nearly quantitatively. 38,146,149,150) The plasmons are generated by on-site charge fluctuations, not by bond-charge ones. The on-site charge excitations are described by the usual density-density correlations and given in the large-N theory by  Figure 12 shows a map of the spectral weight Imχ c (q, q z , ω) in the plane of excitation energy ω and in-plane momentum q along the symmetry axes of (π, π)-(0, 0)-(π, 0)-(π, π).
Below ω ≈ 0.8t, there is a particle-hole continuum coming from individual charge excitations.
π,π) (q x ,q y ) The continuum does not depend much on q z and the result for q z = 0 is presented. We find no strong spectral weight near zero energy, implying that there is no (on-site) charge order tendency, which contrasts with the case of bond-charge orders (see Fig. 11). In a high energy region, there is a sharp peak for q z = 0. This is the well-known optical plasmon and reproduces the early EELS data. 135,136) The plasmon dispersion changes drastically around q = (0, 0) once q z becomes finite. [151][152][153] As a representative, we plot the plasmon dispersion for q z = π in Fig. 12. While the plasmon dispersion remains essentially the same as that for q z = 0 far away from q = (0, 0), the plasmon energy softens substantially near q = (0, 0) and exhibits a strong dispersion there, in sharp contrast to that for q z = 0.
It should be noted that the plasmon energy has a gap at q = (0, 0) for a finite q z . This gap comes from the presence of the finite interlayer hopping integral t z (Ref. 38,133,151,153,154), which was missed in many theoretical studies. [126][127][128][129][130][131][132]152) The magnitude of the gap is proportional to t z in a small t z region and vanishes at t z = 0. In this sense, we call the plasmons for a finite q z as acoustic-like plasmons. The t z dependence of the optical plasmon is almost negligible. 38) A characteristic feature of plasmons lies in the q z dependence of the plasmon energy. Figure 13 shows a map of the spectral weight of plasmons in the plane of q z and ω at a small q. The plasmon energy rapidly decreases with increasing q z and stays almost constant in q z > π/3; this rapid change is pronounced more when a smaller q is chosen. The plasmon intensity, on the other hand, increases with increasing q z , following almost a q 2 z dependence at small q z .
What happens if the long-range Coulomb interaction is replaced by a short-range one, which is more conventional in research of cuprates? In this case, instead of plasmons, the zero-sound mode is realized and its dispersion [ Fig. 14(a)] becomes similar to Fig. 12 around q = (0, 0). A crucial difference appears in the q z dependence. As shown in Fig. 14(b), the zerosound mode energy increases with increasing q z at a small q, which is qualitatively different from the plasmon case shown in Fig. 13. This is because the zero-sound mode becomes gapless at q = (0, 0) and q z = 0 [inset in Fig. 14(a)], whereas the plasmon has a large gap as the optical plasmon (Fig. 12).
It was shown 149) that the acoustic-like plasmon dispersion obtained theoretically well reproduces the peak position of the charge excitations observed in electron-doped cuprates [137][138][139] with doping δ = 0.15, and for hole-doped cuprates 140) with δ = 0.125 and 0.25. Moreover, as seen in Fig. 13, a crucial feature of plasmons is recognized in their characteristic q z dependence. This feature was confirmed not only in electron-doped cuprates 144) but also in hole-doped cuprates. 146) Hence the charge excitations around q = (0, 0) can be summarized as follows. They originate from the acoustic-like plasmons for a finite q z and become the usual optical plasmon for q z = 0. The former explains many RIXS data [137][138][139][140][141][144][145][146][147] and the latter reproduces the early data of plasmons. 135,136) The acoustic-like plasmons have a gap at q = (0, 0), which has not been confirmed in experiments. 146)

Summary and outlook
Given that high-temperature cuprate superconductors are realized by carrier doping into a Mott insulator, the understanding of the charge dynamics is definitely indispensable to the cuprate physics. In this article, we have focused on the electronic nematic order, bond-charge orders, and plasmons.
A crucial insight obtained in theory is that cuprates can be close to the electronic nematic instability. 42,43,155) The nematic order is driven by the J-term and the Coulomb repulsion on a square lattice and tends to be enhanced in the underdoped region 43,[53][54][55] as well as around van Hove filling 43,44) in the hole-doped case; the nematic tendency is weak in the electron-doped case. 61) The nematic order is one of competing orders and can be preempted by antiferromagnetism and superconductivity. 43,58,[62][63][64][65] Nonetheless, even in such a case, the system can retain a large susceptibility to an external xy anisotropy caused by, for example, crystal anisotropy, uniaxial pressure, and external strain. 43) As a result, the system can exhibit a big xy anisotropy in spite of a small external anisotropy. 42, 43, 53-55, 57, 69-72) This physics traces back to a theoretical proposal in 2000 42,43) and is now frequently discussed as evidence of the underlying nematic order. 21,[74][75][76][77][78][79][80]82) The nematic fluctuations generate a pseudogap in the one-particle spectral function in a perturbative calculation to first order. 85) However, more elaborate calculations found no pseudogap. 85) Even in this case, the nematic fluctuations yield a large momentum dependence of the spectral weight with d-wave symmetry along the Fermi surface, leading to a Fermi-arc-like feature. 53,55,85) Given that actual materials frequently contain a small xy anisotropy, the phase diagram shown in Fig. 8 may serve for a basis to discuss a global understanding. 101) Theoretical predictions about the nematic physics are supported in various experiments such as inelastic neutron scattering, 21,[74][75][76] ARPES, 80) Compton scattering, 82) Raman scattering, 81) measurements of Nernst coefficients 77,78) and magnetic torque 79) in hole-doped cuprates. Nonetheless it remains to be studied how the nematic instability and nematic fluctuations can be connected with the pseudogap. The nematic physics is also discussed in ironbased superconductors. 156) The origin of the nematicity, however, may not lie in the dPI [42][43][44] nor charge stripes, 26) but the orbital [157][158][159][160][161] or spin nematicity. 96,97,162,163) This topic was not covered in this article.
Bond-charge orders range from the so-called flux phase 105,106) to s-wave, d-wave, and unidirectional orders. 37,61) These different charge-order tendencies are driven by the spin-spin interaction such as the J-term and can be handled on an equal footing in a large-N theory of the t-J model on a square lattice. In the hole-doped region, a large number of theoretical studies 119,[164][165][166][167][168][169][170][171][172] including the large-N theory 37) were performed, but the consensus has not been obtained on the understanding of the RIXS data. [15][16][17] Given the fact that the charge ordering tendency was observed inside the pseudogap phase, 6) but calculations did not handle the pseudogap appropriately, the pseudogap seems to play an important role to stabilize the charge orders. It is an open issue to pin down the mechanism of the charge order in hole-doped cuprates. On the other hand, in the electron-doped case, where the effect of the pseudogap is weak or absent, 173) the large-N theory can capture the RXS 18,110) and RIXS 111) data very well.
It is d-wave bond-charge order that explains the charge peak observed along the (1, 0) direction for various doping rates. 123,124) The large-N theory also predicts large spectral weight of various bond-charge orders around q = (π, π) (Refs. 122,123), which has not been tested in RXS nor RIXS.
Plasmons originate from the usual on-site charge excitations in the presence of the longrange Coulomb interaction, 125) not from the bond-charge excitations. Moreover the layered structure common to high-T c cuprate superconductors is needed to be considered 38) beyond the widely studied two-dimensional models on a square lattice. The well-known optical plasmon 135,136) is obtained for q z = 0 and the acoustic-like plasmons recently re-ported by RIXS [144][145][146][147] correspond to a finite q z with a V-shape dispersion around q = (0, 0) (Refs. 38,122,149,150). A crucial aspect is the characteristic q z dependence, which serves to identify the origin of the charge excitations. 149) The plasmon scenario can explain the charge excitation spectra observed by RIXS for both hole-(Refs. 146,149) and electron-doped 149,150) cuprates almost quantitatively. Plasmons are the collective on-site charge excitation modes and are realized above the continuum. The continuum spectrum, on the other hand, does not exhibit a strong peak structure at a certain momentum and there is no tendency of the usual charge-density-wave instability. 38) While the subjects we have discussed are limited, we hope that the present article serves for a sound basis toward further experimental and theoretical studies on the origin of the pseudogap and ultimately the high-T c mechanism.

Acknowledgment
The theoretical insights into the electronic nematic order are owed to collaboration with P.