Introduction

Superconductivity is driven by forming Cooper pairs of electrons1. This is achieved even above the boiling point of liquid nitrogen at ambient pressure in cuprate superconductors2. The mechanism of such a remarkable phenomenon has been a central issue in condensed matter physics since their discovery3. The one-particle property of electrons may possess an important hint to solve it. In fact, electrons are not independent inside a material, but acquire the self-energy through various interactions—it is such electrons which drive the superconductivity.

A well-known interaction is electron–phonon coupling, which renormalizes the electron band to yield a kink at the corresponding phonon energy4. A kink was actually observed in not only conventional metals5,6 but also high-temperature cuprates superconductors7,8. Given that the electron–phonon coupling is the conventional mechanism of superconductivity9, its role in the high-Tc mechanism drew much attention. On the other hand, it is also well recognized that the formation of a kink is not a special feature of electron–phonon coupling, but rather a manifestation of coupling to bosonic fluctuations from a general point of view10. Magnetic fluctuations are bosonic fluctuations and can in fact yield a kink in the electron dispersion at the energy of the magnetic resonance mode11,12,13,14.

How about charge fluctuations, which are also bosonic ones? While the understanding of the charge fluctuations is crucial to the high-Tc mechanism in cuprates, it is only recently when the charge dynamics was revealed in momentum-energy space by an advanced technique of resonance inelastic X-ray scattering (RIXS). In particular, low-energy collective charge excitations were revealed in electron-doped15,16,17 and hole-doped cuprates18,19. The characteristic in-plane and out-of-plane dependence allowed to identify these excitations as low-energy plasmons20,21,22,23, which were also discussed for layered metallic systems in the 1970s24,25,26. Dispersing charge modes were reported previously by other groups, too, but interpreted differently, not as plasmons27,28,29,30,31. While the dispersion was presumed to be likely an acoustic mode15,16,18,19, it was found very recently that the low-energy plasmons are gapped at the in-plane zone center for the infinite-layered electron-doped cuprate Sr0.9La0.1CuO217, in agreement with a theoretical study—the gap is predicted to be proportional to the interlayer hopping20. These low-energy plasmons may be referred to as acoustic-like plasmons. While the optical plasmon itself was observed already around 1990 by electron energy-loss spectroscopy32,33 and optical spectroscopy34, these recent advances to reveal the charge excitation spectrum in cuprates attract renewed interest and plasmons offer a hot topic in the research of cuprate superconductors.

Electron–plasmon coupling was studied mainly for weakly correlated materials: alkali metals such as Na and Al35, Bi36, graphene37,38, monolayer transition-metal dichalcogenides (Mo, W)(S, Se)239, semiconductors39,40,41, and diamond41. Here the so-called replica bands are known to be generated by coupling to plasmons; they are also referred to as plasmon satellites41,42,43 especially when momentum is integrated. The spectrum of the replica band is usually very broad and has low spectral weight38,39,40,41,44,45,46, making it difficult to confirm it in experiments. Recently, however, the replica band was successfully resolved in various systems—graphene42,47,48, two-dimensional electron systems49,50, and SrIrO3 films51. The presence of the replica band implies that the one-particle Green’s function has poles. That is, the replica band corresponds to the dispersion relation of quasiparticles dubbed as plasmarons52,53,54,55.

Can we expect plasmarons in high-temperature cuprate superconductors? This is a far from obvious issue. First of all, cuprates are strongly correlated electron systems and thus it is reasonable to distinguish cuprates from weakly correlated systems as the ones discussed in the previous paragraph—a direct analogy between them is not trivial. Second, cuprates are layered systems, where not only the conventional optical plasmon, but also many acoustic-like plasmons are present20. This is also a situation different from previous studies of plasmarons35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51.

In this paper, we show that instead of yielding a kink, plasmons in cuprates lead to plasmarons—similar to weakly correlated systems. A common feature lies in the singularity of the long-range Coulomb interaction in the limit of long wavelength. However, the underlying physics is different. Instead of usual charge density-density correlations, fluctuations associated with the local constraint—non-double occupancy of electrons at any site—are responsible for the emergence of plasmarons. We find that cuprates can host plasmarons near the optical plasmon energy below (above) the quasiparticle band in electron-doped (hole-doped) cuprates.

Results

Analytical scheme

Cuprate superconductors are doped Mott insulators—strong correlations of electrons are believed to be crucial56,57. The t–J model is a microscopic model of cuprates superconductors and is derived from the three-band58 and one-band59 Hubbard model. It reads

$$H=-\mathop{\sum}\limits_{i,j,\sigma }{t}_{ij}{\tilde{c}}_{i\sigma }^{{{{\dagger}}} }{\tilde{c}}_{j\sigma }+\mathop{\sum}\limits_{\langle i,j\rangle }{J}_{ij}\left({\overrightarrow{S}}_{i}\cdot {\overrightarrow{S}}_{j}-\frac{1}{4}{n}_{i}{n}_{j}\right)+\frac{1}{2}\mathop{\sum}\limits_{i\ne j}{V}_{ij}{n}_{i}{n}_{j}\,,$$
(1)

where \({\tilde{c}}_{i\sigma }^{{{{\dagger}}} }\) (\({\tilde{c}}_{i\sigma }\)) are the creation (annihilation) operators of electrons with spin σ(=↑, ↓) in the Fock space without double occupancy at any site—strong correlation effects, \({n}_{i}={\sum }_{\sigma }{\tilde{c}}_{i\sigma }^{{{{\dagger}}} }{\tilde{c}}_{i\sigma }\) is the electron density operator, and \({\overrightarrow{S}}_{i}\) is the spin operator. While cuprates are frequently modeled on a square lattice, we take the layered structure of cuprates into account and consider a three-dimensional lattice to describe plasmons correctly. The sites i and j run over such a three-dimensional lattice. The hopping tij takes the value t\(({t}^{{\prime} })\) between the first (second) nearest-neighbor sites in the plane and is scaled by tz between the planes. The exchange interaction Jij = J is considered only for the nearest-neighbor sites inside the plane as denoted by 〈i, j〉—the exchange term between the planes is much smaller than J (Thio et al.60). Vij is the long-range Coulomb interaction.

It is highly nontrivial to analyze the strong correlation effects systematically. While a variational approach is powerful61, here we employ a large-N technique in a path integral representation in terms of the Hubbard operators62. In the large-N scheme, the number of spin components is extended from 2 to N and physical quantities are computed by counting the power of 1/N systematically. One of the advantages of this method is that it treats all possible charge excitations on an equal footing63,64. There are two different charge fluctuations: on-site charge fluctuations describing usual charge-density-wave and plasmons, and bond-charge fluctuations describing charge-density-waves with an internal structure such as d-wave and s-wave symmetry, including the flux phase. Explicit calculations clarified that those two fluctuations are essentially decoupled to each other65. Since we are interested in plasmons, we focus on the former fluctuations.

Because of the local constraint that double occupancy of electrons is prohibited at any lattice site, the charge fluctuations are described by a 2 × 2 matrix Dab(q, iνn) with a, b = 1, 2; q is the momentum of the charge fluctuations and νn a bosonic Matsubara frequency. While D11 is the usual density-density correlation function, D22 is a special feature of strong correlation effects—it describes fluctuations associated with the local constraint. As we shall clarify, this D22 plays the central role in the formation of plasmarons. Naturally there is also the off-diagonal component D12(=D21).

In the large-N theory, Dab(q, iνn) is renormalized already at leading order. After the analytical continuation iνn → ν + iΓch, where Γch(>0) is infinitesimally small, the full charge excitation spectrum is described by \({{{{{{{\rm{Im}}}}}}}}{D}_{ab}({{{{{{{\bf{q}}}}}}}},\nu )\)—Bejas et al.65 reported a comprehensive analysis of \({{{{{{{\rm{Im}}}}}}}}{D}_{ab}({{{{{{{\bf{q}}}}}}}},\nu )\). In particular, \({{{{{{{\rm{Im}}}}}}}}{D}_{11}({{{{{{{\bf{q}}}}}}}},\nu )\) predicted acoustic-like plasmon excitations with a gap at q = (0, 0, qz) for qz ≠ 0 as well as the well-known optical plasmon at q = (0, 0, 0)20. Close inspections revealed that the predicted plasmon excitations explain semiquantitatively charge excitation spectra reported by RIXS for both hole-doped18,21 and electron-doped17,21,22 cuprates.

Charge fluctuations renormalize the one-particle property of electrons, which can be analyzed by computing the electron self-energy. This requires involved calculations in the large-N theory because one needs to go beyond leading order theory. At order of 1/N, the imaginary part of the self-energy is obtained as66

$${{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )=\mathop{\sum}\limits_{a,b=1,2}{{{{{{{\rm{Im}}}}}}}}{\Sigma }_{ab}({{{{{{{\bf{k}}}}}}}},\omega )\,,$$
(2)

where

$${{{{{\rm{Im}}}}}}{\Sigma }_{ab}({{{{{\bf{k}}}}}},\omega )= \frac{-1}{{N}_{s}{N}_{z}}\mathop{\sum}\limits_{{{{{{\bf{q}}}}}}}{{{{{\rm{Im}}}}}}{D}_{ab}({{{{{\bf{q}}}}}},\nu ){h}_{a}({{{{{\bf{k}}}}}},{{{{{\bf{q}}}}}},\nu ){h}_{b}({{{{{\bf{k}}}}}},{{{{{\bf{q}}}}}},\nu )\\ \times \left[{n}_{{{{{{\rm{F}}}}}}}\left(-{\varepsilon }_{{{{{{\bf{k}}}}}}-{{{{{\bf{q}}}}}}}\right)+{n}_{{{{{{\rm{B}}}}}}}\left(\nu \right)\right].$$
(3)

Here ν = ω − εk−q, εk is the electron dispersion obtained at leading order, ha(k, q, ν) a vertex describing the coupling between electrons and charge excitations, nF and nB the Fermi and Bose distribution functions, respectively, Ns the total number of lattice sites in each layer, and Nz the number of layers; see “Methods” for the explicit forms of Dab(q, ν), εk, and ha(k, q, ν). The real part of Σ(k, ω) is computed by the Kramers-Kronig relations. Since the electron Green’s function G(k, ω) is written as G−1(k, ω) = ω + iΓsf − εk − Σ(k, ω), we obtain the one-particle spectral function \(A({{{{{{{\bf{k}}}}}}}},\omega )=-\frac{1}{\pi }{{{{{{{\rm{Im}}}}}}}}G({{{{{{{\bf{k}}}}}}}},\omega )\):

$$A({{{{{{{\bf{k}}}}}}}},\omega )=-\frac{1}{\pi }\frac{{{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )-{\Gamma }_{{{{{{{{\rm{sf}}}}}}}}}}{{[\omega -{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{{{{{{{\rm{Re}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )]}^{2}+{[{{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )-{\Gamma }_{{{{{{{{\rm{sf}}}}}}}}}]}^{2}}\,,$$
(4)

where Γsf(>0) originates from the analytical continuation in the electron Green’s function.

The quasiparticle dispersion appears as poles of A(k, ω), i.e., a sharp peak structure of A(k, ω), and crosses the Fermi energy. On top of that, A(k, ω) can exhibit other sharp features. If plasmons themselves are responsible for yielding additional poles of A(k, ω), namely fulfilling the condition of \(\omega -{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{{{{{{{\rm{Re}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )=0\) with a relatively small value of \(| {{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )|\), A(k, ω) exhibits a peak describing electrons coupling to plasmons, namely plasmarons52,53,54,55, with a damping controlled by \({{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\). Note that the charge excitation spectrum also contains usual particle-hole excitations, the so-called continuum spectrum, which can also lead to poles in A(k, ω). Thus additional poles of A(k, ω) do not necessarily signal the emergence of plasmarons.

While the t–J model in Eq. (1) contains spin fluctuations, they appear at order of O(1/N) in the present theory whereas charge fluctuations appear at O(1). Hence when we compute the electron self-energy at order of 1/N, only charge fluctuations enter Eq. (3), which is suitable to study the role of plasmons exclusively in the one-particle spectral function.

One-particle spectral function

A choice of model parameters is not crucial to our major conclusions. Here we present results which can be applied directly to electron-doped cuprates, especially La1.825Ce0.175CuO4 (LCCO); details of model parameters are given in “Methods” and Supplementary Note 3.

Figure 1 shows the one-particle spectral function A(k, ω) along the direction (π, π)–(0, 0)–(π, 0)–(π, π); kz dependence is weak and thus kz = π is taken throughout the present paper—the value of kz shall be omitted for simplicity. The spectrum around ω = 0 is the quasiparticle dispersion renormalized by charge fluctuations. In contrast to the case of electron–phonon coupling67,68 and magnetic fluctuations69,70, it does not exhibit a kink structure. Rather it is described by the dispersion εk (solid curve) multiplied by some constant Z(=0.29) as shown by a dotted curve in Fig. 1. This implies that the renormalization factor depends weakly on k and the quasiparticle spectral weight is reduced down to 0.29 by charge fluctuations.

Fig. 1: Intensity map of the one-particle spectral function A(k, ω).
figure 1

The map is focused on a low-energy region along the direction (π, π)–(0, 0)–(π, 0)–(π, π) with kz = π. For comparison, the quasiparticle dispersion obtained at leading order (solid curve) and a renormalized dispersion multiplied by Z = 0.29 (dotted curve) are superimposed. The inset shows A(k,ω) in a wider energy window. The additional band in ω < 0 corresponds to plasmarons, which we shall clarify in the present work.

Full size image

Charge fluctuations also generate additional bands as shown in the inset in Fig. 1. There are two major bands: a low-energy incoherent band near ω ≈ −1.5t and a high-energy side band with a large dispersion in 4 < ω/t < 7—the spectral weight of the former is about 10% and that of the latter is about 60%. The reason to call a side band instead of an incoherent one for the high-energy feature lies in that \({{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) almost vanishes in such a high-energy region, leading to a coherent feature.

Since the low-energy incoherent band disappears when the long-range Coulomb interaction is replaced by a short-range Coulomb interaction—the high-energy one still remains66, a coupling to plasmons is crucial to the low-energy feature. The major point of the present work is to elucidate that the low-energy one corresponds to plasmarons52,53,54,55 and is essentially the same as the so-called replica band discussed in weakly correlated electron systems35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51. In the following we focus on an energy window −2 ≤ ω/t ≤ −1.

Relevant contributions to the formation of plasmarons

As seen in Eq. (2), ImΣ(k, ω) is given by the sum of four components. To elucidate the relevant contribution to forming plasmarons, we may introduce an auxiliary parameter r(≥0) as

$${{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega ;r)={{{{{{{\rm{Im}}}}}}}}{\Sigma }_{11}+r\times \left({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{22}+2\,{{{{{{{\rm{Im}}}}}}}}{\Sigma }_{12}\right)\,,$$
(5)

where the arguments on the right hand side are omitted for simplicity and the fact that ImΣ12 is equal to ImΣ21 was used. The case of r = 1 corresponds to the physical situation seen in Eq. (2). We then compute A(k, ω) for several choices of r in Fig. 2a–d, where a different color scale is used to highlight the weak feature in r < 1; see Supplementary Note 1 for a wider energy window. Upon decreasing r, the incoherent band loses intensity substantially, fades away, and finally becomes invisible in r ≲ 0.4. This clearly indicates that the incoherent band is driven by components involving a, b = 2, namely fluctuations associated with the local constraint—a direct consequence of the strong correlation effect.

Fig. 2: Analysis of plasmarons in terms of Eq. (5).
figure 2

a–d Intensity map of the spectral function A(k, ω) computed with ImΣ(k, ω; r) [Eq. (5)] in −2 ≤ ω/t ≤ −1 for several choices of r along the direction (π, π)–(0, 0)–(π, 0)–(π, π) with kz = π. The plasmaron band in (a) fades away upon decreasing r, indicating that fluctuations associated with the local constraint are crucially important to the plasmarons. Note a different color scale in each panel. e–h Imaginary and real parts of Σ(k, ω) as a function of ω at k = (0, 0) and (π, 0) for several choices of r. The peak of ImΣ(k, ω) in (e) and (g) is determined by the optical plasmon. The line of ω − εk is also shown in (f) and (h). The plasmaron energy is determined by its crossing point of ReΣ(k, ω) on the lower energy side when r is close to 1.

Full size image

The corresponding \({{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) and \({{{{{{{\rm{Re}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) are also shown in Fig. 2e–h for two choices of k. \({{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) exhibits a sharp peak at ω ≈ −1.5t and − 1.3t for k = (0, 0) and (π, 0), respectively, and the peak is suppressed and broadened with decreasing r. The peak structure of \({{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) yields a large dip structure in \({{{{{{{\rm{Re}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) at slightly lower energy than the peak energy of \({{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) via the Kramers-Kronig relations. Consequently, the term \(\omega -{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{{{{{{{\rm{Re}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) in Eq. (4) can vanish at two energies when r is close to 1: one is very close to the peak energy of ImΣ(k, ω) and the other corresponds to the tail of the dip structure of ReΣ(k, ω). Since ImΣ(k, ω) becomes small at the latter energy, A(k, ω) forms a peak there with a damping controlled by ImΣ(k, ω). Hence this peak is incoherent, but is a resonance in the sense that \(\omega -{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{{{{{{{\rm{Re}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )=0\) is fulfilled.

Which one is more crucial to the incoherent band, ImΣ12 or ImΣ22? To answer this, we have studied each component of ImΣab(k, ω) and computed the spectral function for each of them. We can check that ImΣ22 is responsible for the formation of the incoherent band. The component of ImΣ12 works to sharpen the incoherent band by reducing the absolute value of the imaginary part of the self-energy (see Supplementary Note 2 for details).

Role of plasmons

What is then the role of plasmons? Plasmons are described in terms of the charge-charge correlation function, namely poles of \({{{{{{{\rm{Im}}}}}}}}{D}_{11}\) in the present theory. Because of the matrix structure of \({{{{{{{\rm{Im}}}}}}}}{D}_{ab}\), the poles are determined by the zeros of its determinant. Thus all components of \({{{{{{{\rm{Im}}}}}}}}{D}_{ab}\) contain the same poles as those in \({{{{{{{\rm{Im}}}}}}}}{D}_{11}\) and thus describe the same plasmons equally (see Figs. 1 and 8 in Bejas et al.65 for explicit calculations). Because of the layered structure of cuprates, plasmons have various branches depending on the value of qz: the usual optical plasmon corresponds to qz = 0 and the acoustic-like branches to qz ≠ 020. Their energy varies in 0.07 ≤ ν/t ≤ 1.15 around q = (0, 0, qz) for the present parameters. The peak of \({{{{{{{\rm{Im}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )\) at ω = −1.5t (−1.3t) in Fig. 2e [Fig. 2g] is determined by the optical plasmon at ν = 1.15t—the energy difference is easily read off from Eq. (3): the energy of \({{{{{{{\rm{Im}}}}}}}}{D}_{ab}\) is given by ν = ω − εk−q, the optical plasmon is realized at q = 0, ImDab is odd with respect to ν, and thus the peak of ImΣ(k, ω) is shifted by εk−q with q = 0.

The crucial role of the optical plasmon is also confirmed numerically. In Fig. 3, we compute the self-energy ImΣ(k, ω) by removing a region ∣qz∣ ≤ 2π/10 in the qz summation in Eq. (3) so that the contribution from the optical plasmon becomes zero. We then observe that the peak structure in ImΣ(k, ω) completely disappears, not shifts to another energy window. The resulting A(k, ω) no longer forms any structure there.

Fig. 3: Role of the optical plasmon for plasmarons.
figure 3

The self-energy and the spectral function are computed at k = (π, 0) by removing a region ∣qz∣ ≤ 2π/10 in the qz summation in Eq. (3), namely without contributions from the optical plasmon. They do not form any structure. Dotted curves are the corresponding results with the full qz summation. The line of ω − εk is also given. The contrast between the solid and dotted curves demonstrates the importance of the optical plasmon to forming the plasmarons.

Full size image

It might be puzzling—on one hand, the optical plasmon is responsible for the formation of the plasmarons (Fig. 3), but on the other hand, the same plasmons in \({{{{{{{\rm{Im}}}}}}}}{D}_{11}\) and \({{{{{{{\rm{Im}}}}}}}}{D}_{12}\) do not generate them (Fig. 2). The vertex function ha(k, q, ν) in Eq. (3) was checked not to be important. The key lies in the role of the long-range Coulomb interaction Vq, which diverges as q−2 in the limit of q → 0. As is well known, this is the very reason why the plasmons are realized4—the inverse of the charge response function or the determinant of Dab vanishes along the plasmon dispersion. The crucial difference among \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{ab}\) appears in the numerator. To see this we study the explicit form of Dab, which is given by

$${D}_{ab}({{{{{{{\bf{q}}}}}}}},\nu )=\frac{1}{{{{{{{{\mathcal{D}}}}}}}}}\left(\begin{array}{cc}-{{{\Pi }}}_{22}({{{{{{{\bf{q}}}}}}}},\nu )&{{{\Pi }}}_{12}({{{{{{{\bf{q}}}}}}}},\nu )-\frac{N\delta }{2}\\ {{{\Pi }}}_{12}({{{{{{{\bf{q}}}}}}}},\nu )-\frac{N\delta }{2}\quad &-{{{\Pi }}}_{11}({{{{{{{\bf{q}}}}}}}},\nu )+\frac{N{\delta }^{2}}{2}\left({V}_{{{{{{{{\bf{q}}}}}}}}}-{J}_{{{{{{{{\bf{q}}}}}}}}}\right)\end{array}\right)\,.$$
(6)

Here \({{{{{{{\mathcal{D}}}}}}}}\) is the determinant of the matrix \({[{D}_{ab}({{{{{{{\bf{q}}}}}}}},\nu )]}^{-1}\), δ the doping rate, Jq the superexchange interaction in momentum space, Πab a bubble describing particle-hole excitations with appropriate vertex functions ha(k, q, ν), and N the number of spin components (N = 2 corresponds to the physical situations); a complete expression of each quantity is given in “Methods”. In the limit of q → 0, we can obtain by virtue of Vq ~ q−2

$${{{{{{{\rm{Im}}}}}}}}{D}_{22}({{{{{{{\bf{q}}}}}}}},\nu ) \sim \frac{{V}_{{{{{{{{\bf{q}}}}}}}}}^{2}{{{{{{{\rm{Im}}}}}}}}{{{\Pi }}}_{22}}{{({{{{{{{\rm{Re}}}}}}}}{{{{{{{\mathcal{D}}}}}}}})}^{2}+{({{{{{{{\rm{Im}}}}}}}}{{{{{{{\mathcal{D}}}}}}}})}^{2}}\,.$$
(7)

The other components of \({{{{{{{\rm{Im}}}}}}}}{D}_{11}\) and \({{{{{{{\rm{Im}}}}}}}}{D}_{12}\) become smaller by order of \({V}_{{{{{{{{\bf{q}}}}}}}}}^{-2}\) and \({V}_{{{{{{{{\bf{q}}}}}}}}}^{-1}\), respectively. Hence, \({{{{{{{\rm{Im}}}}}}}}{D}_{22}\) becomes dominant over the other components in the limit of q → 0 and thus \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{22}\) has a sizable contribution compared with the other components. This explains the reason why \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{22}\) is responsible for the plasmarons, although all components of \({{{{{{{\rm{Im}}}}}}}}{D}_{ab}\) equally describe the same plasmons.

Discussion

From Eqs. (3) and (7), one may recognize that the mathematical structure of \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{22}\) in the limit of q → 0 is the same as the well-known expression of the self-energy in weak coupling theory,

$${{{{{{{\rm{Im}}}}}}}}{\Sigma }^{{{{{{{{\rm{RPA}}}}}}}}}({{{{{{{\bf{k}}}}}}}},\omega )=\frac{-1}{{N}_{z}{N}_{s}}\mathop{\sum}\limits_{{{{{{{{\bf{q}}}}}}}}}{{{{{{{\rm{Im}}}}}}}}{D}^{{{{{{{{\rm{RPA}}}}}}}}}({{{{{{{\bf{q}}}}}}}},\nu )\left[{n}_{{{{{{{{\rm{F}}}}}}}}}\left(-{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}}^{{{{{{{{\rm{RPA}}}}}}}}}\right)+{n}_{{{{{{{{\rm{B}}}}}}}}}(\nu )\right]\,,$$
(8)

where \({{{{{{{\rm{Im}}}}}}}}{D}^{{{{{{{{\rm{RPA}}}}}}}}}({{{{{{{\bf{q}}}}}}}},\nu )={V}_{{{{{{{{\bf{q}}}}}}}}}^{2}{{{{{{{\rm{Im}}}}}}}}{{{\Pi }}}^{{{{{{{{\rm{RPA}}}}}}}}}({{{{{{{\bf{q}}}}}}}},\nu )\) is the imaginary part of the screened Coulomb interaction computed in the random phase approximation (RPA); \(\nu =\omega -{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}}^{{{{{{{{\rm{RPA}}}}}}}}}\). Therefore weakly correlated electron systems in general can also host plasmarons in principle35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51. However, plasmarons are overdamped in many cases and leave faint spectral weight38,39,40,41,44,45,46. This unfavorable situation is soften when the system has a relatively small band width so that the correlation effect becomes relatively large. In fact, the importance of the small band width to plasmarons was discussed in SrIrO351. See Supplementary Note 5 for explicit results.

In the present t–J model, the band width is very small at order of tδ/2 and δ ≈ 0.1–0.2. Furthermore, the magnitude of all the components of \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{ab}\) is comparable to each other, but the sign of \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{12}\) is the opposite to those of \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{11}\) and \({{{{{{{\rm{Im}}}}}}}}{\Sigma }_{22}\) in ω < 0 (see Supplementary Note 2). Hence after the summation in Eq. (2), \({{{{{{{\rm{Im}}}}}}}}\Sigma\) is substantially reduced in ω < 0. Nonetheless, a small band width allows to fulfill the resonance condition \(\omega -{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{{{{{{{\rm{Re}}}}}}}}\Sigma ({{{{{{{\bf{k}}}}}}}},\omega )=0\) as shown in Fig. 2f and h for r = 1. These features work constructively to host plasmarons in a strongly correlated electron system more than a weakly correlated one.

The dispersion of plasmarons exhibits a dispersive feature similar to the quasiparticle dispersion as shown in Fig. 4—the plasmaron dispersion follows 0.98εk − 1.33t. The value of −1.33t is related to, but not exactly equal to, the optical plasmon energy ν = 1.15t and the factor 0.98 is a renormalization. The plasmaron dispersion can also be fitted to εk − 1.33t approximately. This feature is easily understood intuitively. The energy of the one-particle excitation ω is related to the charge fluctuation energy ν via ω = ν + εk−q in Eq. (3). Since it is the optical plasmon which generates the plasmarons, ν is estimated by its energy and q may be put to zero; recall that ImDab(q, ν) is an odd function with respect to ν. Consequently, the dispersion of plasmarons essentially follows the (bare) quasiparticle dispersion εk. In this sense a term of “replica band” used in weakly correlated electron systems can be inherited even in strongly correlated electron systems.

Fig. 4: Dispersion of plasmarons.
figure 4

It follows 0.98εk − 1.33t (dashed curve) and in this sense, it is a replica band of εk.

Full size image

Figures 1, 2a, and 4 can be applied directly to electron-doped cuprates, especially LCCO. The energy of plasmarons is controlled by the optical plasmon energy, which can be determined precisely by electron energy-loss spectroscopy and optical spectroscopy. Given that the typical energy scale of the optical plasmon in cuprates is around 1 eV, the plasmarons can be tested by angle-resolved photoemission spectroscopy (ARPES) by searching the energy region typically around 1 eV below the electron dispersion especially along the direction (0, 0)–(π, 0) [see the inset of Fig. 1 and Fig. 4]. This energy window has not been studied in detail in ARPES71. Recalling that the plasmarons were detected even in weakly correlated systems such as graphene42,47,48, two-dimensional electron systems49,50, and SrIrO3 films51, there seems a good chance to reveal them also in cuprates.

In experiments72, the optical plasmon energy increases with carrier doping up to 20% doping. Therefore, the energy of plasmarons follows the same tendency. This feature can also be utilized to confirm plasmarons in cuprates. While ARPES is an ideal tool to test plasmarons, X-ray photoemission spectroscopy43 and tunneling spectroscopy42 can also be exploited to detect plasmarons as an emergent satellite peak.

What happens for hole-doped cuprates, where the optical plasmon32,33,34 as well as acoustic-like plasmons17,18,19 typical to layered materials were actually observed? Performing the same analysis as that for the electron-doped cuprates, we predict plasmarons also in hole-doped cuprates. In contrast to the electron-doped case, however, they are realized along the direction (π, 0)–(π, π)–(π/2, π/2) in ω > 0, requiring inverse ARPES to test the plasmaron dispersion; See Supplementary Note 4 for details.

For cuprates, it has been discussed that coupling to bosonic fluctuations yields a kink in the electron dispersion. While plasmons are also bosonic fluctuations, their role in cuprates should be sharply distinguished from phonons7,8 and magnetic fluctuations11,12,13,14. Plasmons do not yield a kink (see Fig. 1), but instead generate plasmarons as an emergent incoherent band (Fig. 4).

The present calculations have been performed in a layered t–J model. If one employs a two-dimensional model, the plasmon dispersion in cuprates cannot be captured especially for the optical plasmon. In this sense, the inclusion of the three-dimensionality of the long-range Coulomb interaction is crucially important to discuss plasmarons in cuprates, although we checked that the interlayer hopping integral tz is not relevant to plasmarons.

A replica band is also discussed in the polar electron–phonon coupling mechanism in TiO273,74,75 and the interplay between the electron–phonon and electron–plasmon couplings was studied76. A clear distinction between those two couplings is made by studying the carrier density dependence of the replica band51. In the present study, however, the electron–phonon coupling is irrelevant because phonon energy is limited below 100 meV in cuprates whereas our relevant energy scale is about 1 eV.

Conclusions

The present large-N theory captures the plasmon excitations observed in both electron- and hole-doped high-temperature cuprate superconductors with a good accuracy so that detailed comparisons with experimental data were made17,18,21,22. We have computed the electron self-energy in the same theoretical framework, but by going beyond leading order theory.

Our major point lies in the indication that cuprates can host plasmarons—quasiparticles coupling to plasmons—near the optical plasmon energy below (above) the quasiparticle dispersion in electron-doped (hole-doped) cuprates; plasmons do not yield a kink in the quasiparticle dispersion, in stark contrast to phonons and magnetic fluctuations. Since plasmarons are found clearly close to momentum (π, 0), where the superconducting gap as well as the pseudogap is enhanced, it is very interesting to explore further the role of plasmarons in the formation of the superconducting gap and the pseudogap in cuprate superconductors.

Our second major point lies in elucidating the mechanism of plasmarons: they are driven by the strong correlation effect—fluctuations associated with the local constraint that imposes no double occupancy of electrons at any site. The underlying physics to generate plasmarons in cuprates is thus different from that in weakly correlated electron systems35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51. However, both have a common mathematical structure to yield plasmarons, establishing a general concept of plasmarons in metals. Plasmarons tend to be well-defined for a system with a smaller band width. This condition is usually fulfilled in cuprates because of strong correlations, but also in a weakly correlated system such as SrIrO3 films51.

Methods

We present a minimal description of the large-N theory of the layered t–J model with the long-range Coulomb interaction—a complete formalism is given in Yamase et al.66.

The electron dispersion εk consists of the in-plane dispersion \({\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}^{\parallel }\) and the out-of-plane dispersion \({\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}^{\perp }\),

$${\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}={\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}^{\parallel }+{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}^{\perp }\,.$$
(9)

At leading order, they are calculated as

$${\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}^{\parallel }=-2\left(t\frac{\delta }{2}+\Delta \right)\left(\cos {k}_{x}+\cos {k}_{y}\right)-4{t}^{{\prime} }\frac{\delta }{2}\cos {k}_{x}\cos {k}_{y}-\mu \,,$$
(10)
$${\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}^{\perp }=-2{t}_{z}\frac{\delta }{2}{\left(\cos {k}_{x}-\cos {k}_{y}\right)}^{2}\cos {k}_{z}\,,$$
(11)

where Δ is the mean value of the bond field, δ the doping rate, and μ the chemical potential. For a given δ, Δ and μ are determined self-consistently by solving the following coupled equations:

$$\Delta =\frac{J}{4{N}_{s}{N}_{z}}\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}}\left(\cos {k}_{x}+\cos {k}_{y}\right){n}_{{{{{{{{\rm{F}}}}}}}}}\left({\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}\right)\,,$$
(12)
$$(1-\delta )=\frac{2}{{N}_{s}{N}_{z}}\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}}{n}_{{{{{{{{\rm{F}}}}}}}}}({\varepsilon }_{{{{{{{{\bf{k}}}}}}}}})\,.$$
(13)

As already mentioned in the Analytical Scheme subsection, charge fluctuations in the t–J model are composed of on-site charge and bond-charge fluctuations. They are, however, essentially decoupled to each other65. Since the former is relevant to the present work, we focus on that. In this case, the bosonic propagator of charge fluctuations, namely Dab(q, iνn), is described by a 2 × 2 matrix with a, b = 1, 2:

$${\left[{D}_{ab}\left({{{{{{{\bf{q}}}}}}}},{{{{{{{\rm{i}}}}}}}}{\nu }_{n}\right)\right]}^{-1}={\left[{D}_{ab}^{(0)}\left({{{{{{{\bf{q}}}}}}}},{{{{{{{\rm{i}}}}}}}}{\nu }_{n}\right)\right]}^{-1}-{{{\Pi }}}_{ab}\left({{{{{{{\bf{q}}}}}}}},{{{{{{{\rm{i}}}}}}}}{\nu }_{n}\right)\,,$$
(14)

where \({D}_{ab}^{(0)}({{{{{{{\bf{q}}}}}}}},{{{{{{{\rm{i}}}}}}}}{\nu }_{n})\) is the bare bosonic propagator,

$${\left[{D}_{ab}^{(0)}\left({{{{{{{\bf{q}}}}}}}},{{{{{{{\rm{i}}}}}}}}{\nu }_{n}\right)\right]}^{-1}=N\left(\begin{array}{cc}\frac{{\delta }^{2}}{2}\left({V}_{{{{{{{{\bf{q}}}}}}}}}-{J}_{{{{{{{{\bf{q}}}}}}}}}\right)&\frac{\delta }{2}\\ \frac{\delta }{2}&0\end{array}\right)\,.$$
(15)

Here \({J}_{{{{{{{{\bf{q}}}}}}}}}=\frac{J}{2}(\cos {q}_{x}+\cos {q}_{y})\) is the superexchange interaction in momentum space and Vq is the long-range Coulomb interaction for a layered system77:

$${V}_{{{{{{{{\bf{q}}}}}}}}}=\frac{{V}_{{{{{{{{\rm{c}}}}}}}}}}{A({q}_{x},{q}_{y})-\cos {q}_{z}}\,,$$
(16)

where \({V}_{{{{{{{{\rm{c}}}}}}}}}={e}^{2}d{(2{\epsilon }_{\perp }{a}^{2})}^{-1}\) and \(A({q}_{x},{q}_{y})=\alpha (2-\cos {q}_{x}-\cos {q}_{y})+1\); e is the electric charge of electrons, a the unit length of the square lattice, d the distance between the layers, α describes the anisotropy between the in-plane and out-of-plane interaction and is given by \(\alpha =\frac{\tilde{\epsilon }}{{(a/d)}^{2}}\) with \(\tilde{\epsilon }={\epsilon }_{\parallel }/{\epsilon }_{\perp }\), where ϵ∥ and ϵ⊥ are the dielectric constants parallel and perpendicular to the planes, respectively. The 2 × 2 matrix Πab is the bosonic self-energy at leading order

$${{{\Pi }}}_{ab}({{{{{{{\bf{q}}}}}}}},{{{{{{{\rm{i}}}}}}}}{\nu }_{n})= -\frac{N}{{N}_{s}{N}_{z}}\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}}{h}_{a}\left({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{q}}}}}}}},{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}}\right)\frac{{n}_{{{{{{{{\rm{F}}}}}}}}}\left({\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}}\right)-{n}_{{{{{{{{\rm{F}}}}}}}}}\left({\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}\right)}{{{{{{{{\rm{i}}}}}}}}{\nu }_{n}-{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}+{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}}}{h}_{b}({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{q}}}}}}}},{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}})\\ -{\delta }_{a1}{\delta }_{b1}\frac{N}{{N}_{s}{N}_{z}}\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}}\frac{{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}-{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}}}{2}{n}_{{{{{{{{\rm{F}}}}}}}}}\left({\varepsilon }_{{{{{{{{\bf{k}}}}}}}}}\right),$$
(17)

and the 2-component vertex is given by

$${h}_{a}({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{q}}}}}}}},\nu )=\left(\frac{2{\varepsilon }_{{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}}+\nu +2\mu }{2}+2\Delta \left[\cos \left({k}_{x}-\frac{{q}_{x}}{2}\right)\cos \left(\frac{{q}_{x}}{2}\right) +\cos \left({k}_{y}-\frac{{q}_{y}}{2}\right)\cos \left(\frac{{q}_{y}}{2}\right)\right],1\right)\,.$$
(18)

We then compute the electron self-energy from charge fluctuations described by Eq. (14) at order of 1/N. This yields Eq. (3) in the main text—its derivation is elaborated in Yamase et al.66.

Fixing temperature to zero, we choose parameters J/t = 0.3, \({t}^{{\prime} }/t=0.3\), tz/t = 0.03, α = 2.9, Vc/t = 18, δ = 0.175, Nz = 10, Γch/t = 0.03, Γsf/t = 0.03, and t/2 = 0.5 eV, which reproduce semiquantitatively the plasmon excitations observed in RIXS for one of the typical electron-doped cuprates LCCO17. The factor of 1/2 in t/2 comes from a large-N formalism where t is scaled by 1/N. We assume N = 2 in comparison with experiments. These parameters were obtained to achieve the best fit to the experimental data under an additional conditions that they should be realistic and do not contradict with the existing knowledge. While Γch and Γsf are positive infinitesimals from the analytical point of view, we employ small, but finite values in actual numerical calculations. This may mimics broadening of the spectrum due to electron correlations at higher orders as well as instrumental resolution. In the figures we presented, all quantities with the dimension of energy are measured in units of t.