Origin of the high-energy charge excitations observed by resonant inelastic x-ray scattering in cuprate superconductors

The recent development of x-ray scattering techniques revealed the charge-excitation spectrum in high-$T_c$ cuprate superconductors. While the presence of a dispersive signal in the high-energy charge-excitation spectrum is well accepted in the electron-doped cuprates, its interpretation and universality are controversial. Since charge fluctuations are observed ubiquitously in cuprate superconductors, the understanding of its origin is a pivotal issue. Here, we employ the layered $t$-$J$ model with the long-range Coulomb interaction and show that an acoustic-like plasmon mode with a gap at in-plane momentum (0,0) captures the major features of the high-energy charge excitations. The high-energy charge excitations, therefore, should be a universal feature in cuprate superconductors and are expected also in the hole-doped cuprates. Acoustic-like plasmons in cuprates have not been recognized yet in experiments. We propose several experimental tests to distinguish different interpretations of the high-energy charge excitations.

to distinguish different scenarios.
We compute the imaginary part of the usual charge susceptibility Imχ c (q, ω) in the layered t-J model with the long-range Coulomb interaction in a large-N scheme (see Methods and Supplemental Material). We show in Fig. 1 the high-energy collective dispersions, namely the peak positions of Imχ c (q, ω), for e-cuprates with doping x = 0.15 (black line), and for h-cuprates with x = 0.125 (green line) and x = 0.25 (purple line). These excitations correspond to plasmons realized in the layered system with a gap at q = (0, 0). For comparison, we include in Fig. 1 the peak position of the charge excitations obtained in the experiments [18][19][20][21] . The agreement with the experimental data is very good in |q | 0.5π for both e-and h-cuprates.
Our obtained dispersion for x = 0.25 has higher energy than that for x = 0.125. This feature captures the experimental data shown in Fig. 1. In fact, the analysis in Ref. 21 finds that the energy at q = (0.46π, 0) increases by a factor of 1.16 when doping is increased from x = 0.125 to x = 0.25. In the present theory, we obtain a factor of 1.30. Similar results were also obtained in the density-matrix renormalization-group calculations in the three-band Hubbard model 21 .
For x = 0.15 in e-cuprates in the large q region, i.e., (0.5π, 0)-(π, 0), the experimental data deviate downward from the dispersion proposed in Ref. 18 (dashed line in Fig. 1) and tend to be closer to our obtained dispersion. Still, the deviation between the experimental data and our results seems substantial, compared with the agreement in the small q region.
However, we think that such a deviation could be related to, as we shall discuss later, the broad spectrum observed in the experiments especially in a large q region.
In addition to the agreement with the experimental data in Fig. 1, the present theory implies the following: i) The high-energy charge excitations correspond to a plasmon mode with a gap at q = (0, 0); the gap is proportional to t z (Ref. 24). ii) Our high-energy charge excitations are present in both e-and h-cuprates. iii) The dispersion around q = (0, 0) has a larger slope in e-cuprates than h-cuprates, consistent with the observation in Ref. 21. iv) The dispersion is rather symmetric between the direction (0, 0)-(π, 0) and (0, 0)-(π, π).
Since the plasmon is a collective mode, it can form a very sharp peak in q -ω space as shown in Fig. 2a. On the other hand, the experiments of Refs. 19-21 do not show a peak signal at q = (0, 0) and in addition, the spectrum is broad and becomes broader with increasing q toward the BZ boundary (Fig. 2d). These features are not seen in Fig. 2a because it was computed in an ideal situation by taking Γ = 0.001t, i.e., the damping is assumed to be very small (see Methods for the definition of Γ). Considering a realistic situation, we compute the charge-excitation spectrum by employing a large Γ. In principle, Γ would depend on momentum and energy, but we take a constant Γ = 0.7t as the simplest case. As shown in Fig. 2c, the spectrum is substantially broadened and becomes broader with increasing q . In addition, the spectrum near q = (0, 0) becomes poorly resolved.
These features are very similar to the experimental results (Fig. 2d). In Fig. 2d, there is strong intensity around q = (0, 0) and ω = 2 eV, which comes from the charge transfer excitations between oxygens and coppers, namely interband excitations. This feature is beyond the scope of the analysis of the present one-band model.
The inclusion of a large Γ is actually invoked theoretically when the spectral line shape is compared with experiments 21 . In Ref. 24 a finite Γ was also used to discuss the temperature dependence of spectral weight 18 . Physically there should be two different broadenings, intrinsic and extrinsic ones. The extrinsic broadening is due to the instrumental resolution, which is about 250 meV (Γ ∼ 0.35t) in Ref. 19, and 130 meV (Γ ∼ 0.2t) in Ref. 18. Because Γ = 10 -3 t (π,π) (0,0) (π,0) (q x ,q y ) We have also calculated the peak area at a given q along the (π, π)-(0, 0)-(0, π) direction and compare it with the experimental results in Fig. 2b. This agreement with the experiment strengthens the idea that the high-energy charge excitations are plasmons.
Our plasmon mode should not be confused with usual optical plasmons, which are actually observed in optical measurements 27 and electron energy-loss spectroscopy 28,29 in cuprate superconductors. The optical plasmon mode is in fact reproduced in our theory by invoking q z = 0 (Ref. 24). However, q z is usually finite in RIXS. Once q z becomes finite, the optical plasmon energy is substantially suppressed to be proportional to the interlayer hopping t z , yielding acoustic-like plasmons as shown in Fig. 1. While this strong q z dependence was, in part, already discussed in Ref. 24, as well as in early theoretical works where t z = 0 was assumed in a layered model 30-32 , we present further results. Figure 3 shows a map of the spectral weight of plasmons in the plane of q z and ω for several choices of Γ 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 more pronounced when a smaller q is chosen. The plasmon intensity, on the other hand, increases with increasing q z , following nearly a q 2 z dependence at small q z . Those qualitative features are independent of the broadening Γ. However, the peak intensity at a small q is suppressed substantially with increasing Γ (see also Fig. 2c).
Hence the q z dependence of plasmons may be well observed for a small Γ. Although the importance of the q z dependence of plasmons was not recognized in experimental papers 18 When q z = 0, V (q) is singular at q = (0, 0), which leads to usual optical plasmons.
However, due to the anisotropy between q and q z , the plasmon energy becomes different when q z is reduced to zero at q = (0, 0). In particular, the plasmon energy would become zero if the interlayer hopping is neglected. This is the reason why the plasmon energy becomes sensitive to the value of q z , especially in a region of a small q .
We have shown that acoustic-like plasmon excitations with a gap at q = (0, 0) due to a finite interlayer hopping describe the main features observed by different experimental for details). While the peak position depends on choices of Γ, its q z dependence is almost negligible. This feature is qualitatively different from plasmons shown in Fig. 3 and thus serves to clarify the underlying physics of the high-energy charge excitations.
One may wish to consider a scenario without the long-range Coulomb interaction, which may replace plasmons with a zero-sound mode in the t-J model. To demonstrate this, we have computed charge excitations in our model by using the short-range Coulomb interaction instead of the long-range one (see Supplemental Material). Our obtained spectrum is shown in Fig. 5a, which is qualitatively similar to Fig. 2a. While it is clear theoretically that the zero-sound mode is fundamentally different from plasmons 34 , their distinction is less clear from an experimental point of view. Hence, we have computed the q z dependence of the zero-sound mode for a small q in Fig. 5b. The zero-sound energy increases with increasing q z in a small q z region, which is qualitatively different from the plasmon case shown in Fig. 3. This is because the zero-sound mode becomes gapless at q = (0, 0) and q z = 0 as shown in the inset of Fig. 5a. Therefore, besides Ref. 33, additional experimental data about the q z dependence of the high-energy charge excitations may confirm the importance of the long-range Coulomb interaction in the charge dynamics in cuprates.
We have demonstrated that acoustic-like plasmon excitations can consistently explain experimental data obtained by different groups [18][19][20][21] . Conceptually plasmons are well known in solids, but the presence of the acoustic-like plasmon mode in cuprates has not been recog-nized yet in experiments. Thus our theoretical recognition of the acoustic-like plasmons in cuprates highlights the importance of charge dynamics in cuprates. Recalling that cuprates have been studied largely by focusing on spin degrees of freedom, it is worth exploring unresolved issues in cuprates such as the origin of the pseudogap and the mechanism of high-T c superconductivity in terms of charge degrees of freedom, including the present acoustic-like plasmons.

Methods
We employ the t-J model on a square lattice by including the interlayer hopping and the long-range Coulomb interaction: wherec † iσ andc iσ are the creation and annihilation operators of electrons with spin σ in the Fock space without double occupancy, and i and j run over a three-dimensional lattice.
iσc iσ is the electron density operator and S i is the spin operator. The hopping t ij takes a value t (t ′ ) between the first (second) nearest-neighbors sites on the square lattice, and t z between the layers. i, j denotes a nearest-neighbor pair of sites. We neglect the magnetic exchange interaction between the planes, which is much smaller than in-plane J (Ref. 35). V ij is the long-range Coulomb interaction on the lattice and it is given in momentum space by 36 where V c = e 2 d(2ǫ ⊥ a 2 ) −1 and where α = (d/a) 2 ǫ /ǫ ⊥ , and ǫ and ǫ ⊥ are the dielectric constants parallel and perpendicular to the planes, respectively; e is the electric charge of electrons; a is the lattice constant in the planes and the in-plane momentum q = (q x , q y ) is measured in units of a −1 ; similarly d is the distance between the planes and the out-of-plane momentum q z is measured in units of d −1 .
Since the Hamiltonian (1)   q z is usually finite. We thus first present results for q z = π as representative ones in Figs. 1 and 2, and then study their q z dependence. The temperature is set to zero.
In the comparison with the experimental data in Fig. 1, we have used t = 750 meV for  In this supplemental material we present i) essential part of our formalism, ii) analysis of individual charge excitations, and iii) short-range Coulomb interaction.

I. THEORETICAL SCHEME
In the path integral formalism 1 , the Hamiltonian (1) can be written in terms of an effective model where fermionic fields interact with the six-component bosonic field Here δR i describes on-site charge fluctuations and is related to i is the Hubbard operator 2 associated with the number of holes at a site i; x is the doped carrier density per site; the factor N comes from the sum over the N fermionic channels after the extension of the spin index σ from 2 to N. δλ i describes fluctuations of the Lagrangian multiplier introduced to impose the constraint of non-double occupancy at any site. r  ab (q, iω n ) of the bosonic field δX a is given by where J(q) = J 2 (cos q x + cos q y ) and the matrix indices a and b run from 1 to 6; q is a three dimensional wavevector and ω n is a bosonic Matsubara frequency.
At leading order, the bare bosonic propagator is renormalized to be where Π ab (q, iω n ) is the 6 × 6 bosonic self-energy due to the coupling between the bosonic field and fermionic fields; n F is the Fermi-Dirac distribution function. The six components interaction vertex is given by The electronic dispersion ǫ k is defined as where the in-plane dispersion ε k and the out-of-plane dispersion ε ⊥ k are given, respectively, by and µ is the chemical potential. Note that the bare hopping integrals t, t ′ , and t z are renormalized by a factor x/2. No incoherent self-energy effects enter the fermionic dispersion at leading order. Note that k z and q z dependences enter only through ǫ k−q in the first column in Eq. (8), whereas the other columns contain only the in-plane momentum q . In Eq. (7), N s and N z are the total number of lattice sites on the square lattice and the number of layers along the z direction, respectively.
All possible charge excitations in the layered t-J model are contained in D ab (q, iω n ) [Eq. (6)] and can be treated on equal footing in the present theoretical scheme 3 . Usual charge fluctuations, namely χ c (r i − r j , τ ) = T τ n i (τ )n j (0) , are associated with the element (1, 1) of the full 6 × 6 D ab in Eq. (6) and is computed in q -ω space as 3,4 χ c (q, iω n ) = −N x 2 2 D 11 (q, iω n ) .
On the other hand, the elements from 3 to 6 of the matrix D ab describe low-energy charge excitations associated with the charge order phenomenon 3 . Since we are interested in the high-energy charge excitations, we focus here on Eq. (12).

II. INDIVIDUAL CHARGE EXCITATIONS
The charge susceptibility χ c (q, iω n ) [Eq. (12) Hubbard X-operators X p0 i and X 0p i by using the two constraints such as X 00 i + p X pp i = N 2 and X pp ′ i = X p0 i X 0p ′ i X 00 i (see Ref. 4). Therefore after the analytical continuation, we have computed ImΠ 22 (q, ω) in Fig. 4. We have checked that the spectrum of ImΠ 22 (q, ω) for a small Γ is indeed similar to the continuum spectrum of Imχ c (q, ω) obtained in Fig. 2a, where the plasmon mode is well separated from the continuum spectrum.

III. SHORT-RANGE COULOMB INTERACTION
As a typical short-range Coulomb interaction, we may take V (q) = V 1 (cos q x + cos q y ) + V 2 cos q z .
The charge excitation spectrum for V 2 = 0 and t z = 0 was already shown in Ref. 5, where a zero-sound mode is realized as collective excitations, which are gapless at q = (0, 0); see