# Fileset

[main.pdf](https://mdr.nims.go.jp/filesets/c2eb631a-2181-40ff-9cbc-0b7287878540/download)

## Creator

[Hiroyuki Yamase](https://orcid.org/0000-0003-0328-5657), [Luciano Zinni](https://orcid.org/0000-0002-9165-9251), [Matías Bejas](https://orcid.org/0000-0003-4254-0622), [Andrés Greco](https://orcid.org/0000-0001-5958-5080)

## Rights

©2026 American Physical Society[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Strong-coupling theory of bilayer plasmon excitations](https://mdr.nims.go.jp/datasets/0013807e-b287-4f54-a1ea-e500f9a11818)

## Fulltext

Strong-coupling theory of bilayer plasmon excitationsHiroyuki Yamase1, Luciano Zinni2, Mat́ıas Bejas3, and Andrés Greco31Research Center of Materials Nanoarchitectonics (MANA),National Institute for Materials Science (NIMS), Tsukuba 305-0047, Japan2Facultad de Ciencias Exactas, Ingenieŕıa y Agrimensura (UNR-CONICET),Avenida Pellegrini 250, 2000 Rosario, Argentina3Facultad de Ciencias Exactas, Ingenieŕıa y Agrimensuraand Instituto de F́ısica Rosario (UNR-CONICET),Avenida Pellegrini 250, 2000 Rosario, Argentina(Dated: October 14, 2025)AbstractRecently plasmon excitations in bilayer lattice systems were studied extensively in the weak-coupling regime. Unlike single-layer systems, these bilayers exhibit two distinct modes, ω±, whichshow characteristic dependences upon the momentum and hopping integrals along the z direction.To apply them to cuprates, strong correlation effects should be considered, but a comprehensiveanalysis has not yet been investigated. In this work, we present a strong-coupling theory to analyzethe charge dynamics of a bilayer system, utilizing the t-J-V model, which includes the long-rangeCoulomb interaction, V , on a lattice. Although our theoretical framework is fundamentally differentfrom the weak-coupling approach, we find that resulting plasmon excitations are similar to thoseof a weak-coupling theory. A key distinction is that our strong-coupling framework reveals anoticeable suppression of particle-hole excitations, which allows the plasmon modes to remainwell-defined over a wider region of momentum. We suggest that the experimentally reportedplasmon excitations in Y-based cuprates can be described by the ω− mode, although we call formore systematic experiments to verify this.1I. INTRODUCTIONThe parent compounds of cuprate superconductors are widely known to be antiferromag-netic Mott insulators. Upon carrier doping, the charge degrees of freedom become active andthe system transitions into a metallic state, which suppresses the antiferromagnetic order. Ahigh-temperature superconducting state emerges at a carrier doping level of approximately5 % for hole-doped and 10–15 % for electron-doped materials, reaching a maximal Tc around16 % doping [1].It is well established that electrons within the CuO2 layers play a central role in high-temperature superconductivity. For this reason, the two-dimensional t-J and Hubbard mod-els are considered the minimal theoretical framework [2]. While the importance of spin-spininteraction is frequently emphasized, the charge degrees of freedom should be equally crucialto understanding the cuprate physics. The full charge dynamics in momentum-energy spacehas recently been revealed comprehensively through the advent of the resonant inelasticx-ray scattering (RIXS) technique [3, 4].In the charge excitation spectrum, two distinct features have been observed. The firstfeature is low-energy excitations around in-plane momentum q∥ = (0.6π, 0) in hole-dopedcuprates [5–21] and around (0.5π, 0) in electron-doped cuprates [22–24]. The second featureis a distinct V-shaped dispersion centered at q∥ = (0, 0). While its origin was initiallydebated [25–30], it is now consistently understood as an acousticlike plasmon, which isa characteristic feature of layered systems [28, 31]. A particular important feature of theplasmon is the rapid decrease of its energy with increasing momentum transfer perpendicularto the layers, qz [31, 32]. Additionally, a gap at q∥ = (0, 0) has been observed [33], which isproportional to the interlayer hopping tz [28]. This dependence provides a valuable way toextract the value of tz, a parameter that is difficult to determine through other experimentaltechniques.The present paper focuses on the plasmon excitations in cuprate superconductors. Whilemost experimental studies have been limited to single- and infinite-layer systems, wherethe unit cell contains one CuO2 plane, a theoretical description requires two additionalfactors beyond the standard t-J and Hubbard models: the long-range Coulomb interaction(LRC)—which in continuum space has a 1/r dependence—and the interlayer hopping, bothof which are essential to correctly capture the observed qz dependence of the plasmon energy2[28, 32, 34].A crucial theoretical question arises regarding how plasmon excitations change with anincreasing number of layers per unit cell. Pioneering studies by Fetter [35] and Griffinand Pindor [36] for a layered electron gas model showed the existence of two modes, oneof which has very low energy. Concurrently, the superconducting onset temperature Tcincreases substantially when the number of CuO2 planes is increased to two or three, up to140 K; for four or more layers, Tc plateaus around 110 K [37].Experimental studies of low-energy plasmon modes for multilayer cuprates remain limited,with a single report in bilayer Y-based compounds [38]—which has been extensively analyzedwithin the weak-coupling random phase approximation (RPA) [38–40]—and a subsequentobservation in trilayer Bi-based cuprates [41].The aim of the present work is to formulate a strong-coupling theory for the bilayer sys-tem. This work complements the very recent weak-coupling analysis in Ref. [39], allowingus to systematically clarify the similarities and differences between the two theoretical ap-proaches. This point is important because some argue that weak-coupling approaches areinadequate for cuprates, given their strongly correlated nature. In fact, alternative frame-works have been proposed, including holographic descriptions of charge dynamics [42–44].Furthermore, the formalism developed here acquires renewed relevance following the recentdiscovery of high-Tc superconductivity in bilayer nickelates [45, 46].The remainder of this paper is organized as follows. In Sec. II, we formulate a large-Ntheory of the layered t-J-V model. The LRC is treated on a lattice, respecting the bilayerstructure [39], rather than a continuum form used in a layered electron gas model [35, 36].In Sec. III, we present our results for charge excitation spectra, which may be comparedto those obtained in the weak-coupling theory [39]. We also investigate the dependenceof these excitation on the LRC, V . Discussions and conclusions are given in Sec. IV andV, respectively. In Appendix A, the complete formalism of our strong coupling theory ispresented. In Appendix B, we discuss whether the collective modes obtained within the t-Jmodel without LRC can account for the experimental data.3II. FORMALISMWe begin with our theoretical analysis by defining the bilayer t-J-V model on a squarelattice. The Hamiltonian is given by:H =−∑i,j,σ,α,βtαβij c̃†iσ,αc̃jσ,β − µ∑i,αni,α + J∑⟨i,j⟩,α(S⃗i,α · S⃗j,α − 14ni,αnj,α)+J⊥2∑i,α̸=β(S⃗i,α · S⃗i,β −14ni,αni,β)+12∑i̸=j,α,βV αβij ni,αnj,β , (1)where i and j run over the three-dimensional lattice sites, α, β = 1, 2 denote the planewithin a unit cell, and c̃†iσ,α (c̃iσ,α) is the creation (annihilation) operator of an electron withspin σ(=↑, ↓) at site i in layer α. The hopping integrals tαβij extend up to third-nearestneighbors on the square lattice of each layer, denoted as t, t′ and t′′, respectively. Alongthe z direction, the hopping within a bilayer is denoted as tz (intrabilayer hopping), andthe hopping between bilayers is given by t′z (interbilayer hopping)—see Fig. 1. J is thestrength of the in-plane spin exchange interaction between nearest-neighbor site ⟨i, j⟩; J⊥is the out-of-plane spin exchange and considered only within the intrabilayer; V αβij is thethree-dimensional LRC. ni,α =∑σ c̃†iσ,αc̃iσ,α is the electron density operator at site i in layerα, S⃗i,α is the spin operator, and µ is the chemical potential. In the t-J-V model, all operatorsare defined in the Fock space without double occupancy, which yields the local constraint:∑σc̃†iσ,αc̃iσ,α ≤ 1 (2)for any site i and layer α.Here we employ a large-N technique in a path integral formalism in terms of the Hubbardoperators [47]. In this scheme, the number of spin component is extended from 2 to N andphysical quantities are computed by the power of 1/N systematically.Setting tz = t′z = J⊥ = V αβij = 0 (α ̸= β) in Eq. (1) reduces the system to two decoupledt-J-V models. Such a model was successfully applied to analyzing the charge dynamics incuprates [28, 32–34, 48–51]. We extend this formalism to the bilayer lattice shown in Fig. 1,which was originally introduced to analyze bilayer nickelates [52]. Our interest lies in bilayerplasmon excitations and we now summarize the key steps of this formalism. The completeformalism is given in Appendix A.At leading order in the bilayer system, the electron Green’s function is obtained as 2× 24FIG. 1. Schematic of the bilayer lattice and corresponding hopping integrals. Each layer forms asquare lattice with lattice constant a; d is the intrabilayer distance and c is the lattice constantalong the z direction.matrix:G(0)αβ(k, iνn) = iνn − ε∥k −ε⊥k eikzdc−ε⊥∗k e−ikzdc iνn − ε∥k−1, (3)where νn is a fermionic Matsubara frequency, and the electron dispersions are obtained asε∥k = −2(tδ2+ χ)(cos kx + cos ky)− 4t′δ2cos kx cos ky − 2t′′δ2(cos 2kx + cos 2ky)− µ , (4)ε⊥k = −[tzδ2(cos kx − cos ky)2 + χ′]− t′zδ2(cos kx − cos ky)2 e−ikz . (5)Here kx, ky and kz are in units of the inverse of the lattice constant a, a, and c, respectively.For a given doping δ, the chemical potential µ and the values χ and χ′ are determinedself-consistently by solving:χ =J8Ns∑k(cos kx + cos ky)[nF (ε1k) + nF (ε2k)], (6)χ′ = − J⊥4Ns∑kε⊥k eikzdc|ε⊥k |[nF (ε1k)− nF (ε2k)], (7)and1− δ =1Ns∑k[nF (ε1k) + nF (ε2k)], (8)where Ns is the number of lattice sites; nF is the Fermi distribution function; from thedeterminant of Eq. (3), the bonding and antibonding bands (α = 1, 2, respectively) can beobtained asεαk = ε∥k − (−1)α|ε⊥k | . (9)5Charge fluctuations are described by the 14 × 14 matrix of the boson propagator atthe order of the 1/N in the bilayer system. But the on-site charge fluctuations includingplasmons that we are interested in are described by a 4× 4 reduced matrix,D−1ab (q, iωn) =[D(0)ab (q, iωn)]−1− Πab(q, iωn) , (10)where a and b run from 1 to 4; q is a three-dimensional vector; ωn is a bosonic Matsubarafrequency. D(0)ab (q, iωn) is the bare bosonic propagator and is obtained as[D(0)ab (q, iωn)]−1= Nδ22[V (q)2− J(q)]δ2δ22[V ′(q)2− J ′(q)]0δ20 0 0δ22[V ′∗(q)2− J ′∗(q)]0 δ22[V (q)2− J(q)]δ20 0 δ20 , (11)where J(q) = (J/2)(cos qx + cos qy) (qx and qy are in units of the inverse of the latticeconstant a) and J ′(q) = (J⊥/4)e−iqzdc (qz is in units of the inverse of the lattice constantc). V (q) and V ′(q) are the intralayer and interlayer Fourier components of the LRC in thebilayer lattice, respectively, which are given byV (q) =VcdetṼ[α (2− cos qx − cos qy)−12h3 −12h1 cos qz], (12)V ′(q) =12VcdetṼ{h2 cos(qzdc)+ h4 cos[qz(1− dc)]− ih2 sin(qzdc)+ih4 sin[qz(1− dc)]}, (13)detṼ = [α (2− cos qx − cos qy)]2 − α (2− cos qx − cos qy) (h1 cos qz + h3)+6c2(c− d)(2c− d)(1− cos qz) . (14)Here Vc =e2c2a2ε⊥, α =c2ε∥a2ε⊥(ε∥ and ε⊥ are the parallel and perpendicular dielectric constants,respectively; e is the electric charge of electron), and h1 =2c(c−2d)(2c−d)(c−d), h2 =2cc−d, h3 = − 4cc−d,and h4 = 2c(c+d)(2c−d)(c−d). Equations (12) and (13) were first derived in Ref. [39] and theyrepresent the only known analytical expressions for the LRC in the bilayer structure. Theself-energy components Πab(q, iωn) in Eq. (10) are calculated in Appendix A.In the large-N formalism for bilayer systems, the charge-charge correlation function χ(ri−rj, τ) = ⟨Tτni(τ)nj(0)⟩ is related in q-ω space to the dressed bosonic propagator Dab(q, iωn)asχ (q, iωn) =N2(δ2)2[D11 (q, iωn) +D33 (q, iωn) +D13 (q, iωn) +D31 (q, iωn)] . (15)6The elements D11 and D33 (D13 and D31) give the intralayer (interlayer) contributions toχ (q, iωn).After performing the analytical continuation iωn → ω+iΓ and adopting the physical valueN = 2, we obtain the imaginary part of the charge-charge correlation function, Imχ(q, ω),which can be directly compared with RIXS measurements. The parameter Γ may accountfor both experimental resolution and spectral broadening from electron correlations [53].III. RESULTSWe present our results organized into three subsection. Focusing on cuprates, we adoptthe following parameter set without losing generality: t′ = −0.3t, t′′ = 0.15t, J = 0.3t,J⊥ = 0 [54]. The lattice parameters are set to a = 3.88 Å, c = 11.68 Å, d = 3.36 Å[38]. Thedielectric constants are chosen as ϵ∥ = 6ϵ0, ϵ⊥ = 2.25ϵ0. These values yield Vc = 38.75t andα = 24. We take hole doping rate δ = 0.21, corresponding to the electron density n = 0.79.The broadening parameter Γ is set to 0.01t except for Fig. 12. We take a temperature ofT = 0.04t to avoid bond-charge instabilities [55, 56]. Our analysis first focuses on the caseof t′z = 0 before examining its effect. Unless stated otherwise, t is the unit of energy.A. Charge excitation spectra with t′z = 01. Overall spectrumFigure 2(a) presents the charge excitations spectrum across a broad range of in-planemomentum q∥ and energy transfer ω. The white dotted curve shows the upper boundaryof the particle-hole continuum. While weak signals are visible below this boundary, themost significant spectral weight is found above it, corresponding to plasmon excitations.These features are qualitatively similar to those observed in a single-layer system (see Fig. 1in Ref. [28]). A crucial distinction emerges at low energy (ω < 0.3t) and small in-planemomentum, specifically around q∥ = (0, 0), where two distinct plasmon modes are presentfor a fixed qz = π. Upon closer inspection, the lower branch is found to be a gaplessplasmon mode that extends into the continuum, a feature characteristic of bilayer systems[39]. Following the nomenclature of Griffin and Pindor [36], we refer to the higher-energyand lower-energy branches as ω+ and ω− modes, respectively.700.51.01.52.0(π,π) (0,0) (π,0) (π,π)ω / tq||<−1.7>0log10|Imχ(q,ω)|(a) qz = π   tz = 0.1t   tz’ = 0qxqyππ00246810(π,π) (0,0) (π,0) (π,π)ω / tq||<−3>0.5log10|Imχ(q,ω)|(b) Weak−coupling RPAFIG. 2. (a) Intensity map of charge excitation spectrum log10 |Imχ(q, ω)| in the plane of in-planemomentum q∥ and energy transfer ω. The white dotted curve represents the upper boundary ofthe particle-hole continuum excitations. The strong intensity corresponds to plasmon modes. (b)The corresponding map in a weak-coupling analysis [39] for the same parameters as (a) except thatVc = 130 is taken.The corresponding weak-coupling RPA result [39] is shown in Fig. 2(b) for exactly thesame parameters as Fig. 2(a) except for a value of Vc so that two plasmon branches are wellvisible. A comparison between Fig. 2(a) and (b) reveals that strong electronic correlationssuppress the particle-hole continuum substantially, allowing the plasmon branches to remainwell defined along the entire Brillouin-zone path, whereas in the weak-coupling schemethey are confined to the vicinity of q∥ = (0, 0) and merge into the continuum away fromq∥ = (0, 0).The remainder of this section will focus on the detailed properties of ω+ and ω− modes.2. Plasmon excitationsTo investigate the plasmon excitations in more detail, we analyze the q-ω map aroundq∥ = (0, 0) for a sequence of qz values as shown in Figs. 3(a)–(e). At qz = 0 [Fig. 3(a)],only the ω+ mode is present. The charge conservation makes the spectral weight vanish atq∥ = (0, 0) because charge fluctuations between the two layers are in-phase. This confirmsthat the ω+ mode corresponds to the well-known optical plasmon. By symmetry, the ω−mode corresponds to out-of-phase fluctuations and is not present for qz = 0, a featurethat has the following physical interpretation: When applying infinitesimally small external800.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(a) qz = 0 tz = 0.1t   tz’ = 000.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(b) qz = 0.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(c) qz = π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(d) qz = 1.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(e) qz = 2πFIG. 3. Intensity maps of log10 |Imχ(q, ω)| for a sequence of qz around a region of q∥ = (0, 0) fortz = 0.1t: (a) qz = 0, (b) qz = 0.5π, (c) qz = π, (d) qz = 1.5π, and (e) qz = 2π. The white dottedcurve denotes the upper boundary of the particle-hole continuum. It goes to zero at q∥ = (0, 0)and qz = 0 in (a).electric field with qz = 0, this external field is strictly uniform along the out-of-plane direction900.20.40.60.80 1 2 3 4 5 6 7 8ω / tqz/π<-2>0log10|Imχ(q,ω)|(a) q|| = (0.02π, 0) tz = 0.1ttz’ = 000.20.40.60.80 1 2 3 4 5 6 7 8ω / tqz/π<-2>0log10|Imχ(q,ω)|(b) q|| = (0.05π, 0)FIG. 4. qz dependence of ω+ mode (higher energy) and ω− mode (lower energy) at (a) q∥ =(0.02π, 0) and (b) q∥ = (0.05π, 0) for tz = 0.1t. The white dotted line is the upper boundary ofthe continuum spectrum and exhibits a sharp drop at qz = 0 because of the vanishing of the ω−mode there.and thus cannot couple to the out-of-phase charge oscillation, i.e., ω− mode.Although the spectral weight at q∥ = (0, 0) vanishes for qz = 0, a finite intrabilayerhopping tz results in a finite energy for the upper boundary of the continuum at q∥ = (0, 0)for qz ̸= 0. This allows for the presence of two plasmon modes, the ω+ and ω− modes, as seenin Figs. 3(b)–(d). These modes should not be confused with the even and odd modes, sincethe LRC couples fluctuations across all layers. Notably, the ω− mode exhibits a gaplessdispersion, which is particularly evident as it enters the continuum around q∥ = (0, 0).However, a significant change occurs at qz = 2nπ where n ̸= 0 [Fig. 3(e)], as the ω−mode becomes gapped. A comparison between Figs. 3(a) and (e) reveals that the plasmonexcitations at qz = 0 are not generic. Instead, the behavior at qz = 2nπ with n ̸= 0 is morerepresentative: both ω± modes are present and gapped, but the intensity of the ω+ modevanishes at q∥ = (0, 0) while it does not for the ω− mode.We next examine the qz dependence of the ω± modes more closely. Figure 4(a) showsresults at a fixed in-plane momentum of q∥ = (0.02π, 0) for tz = 0.1t. The white dotted lineindicates the upper boundary of the particle-hole continuum. The ω+ mode is always presentabove the continuum. Its energy sharply peaks at qz = 2nπ (corresponding to the opticalplasmon) and quickly decreases as qz moves away from these values, remaining largely qzindependent until the next peak. The ω− mode, which is absent at qz = 0, gains intensityimmediately as qz increases. It shows a dispersive feature between ω = 0.06–0.21t with apeak at qz = 2nπ where n ̸= 0. This dispersion is particularly interesting since it occurs by10crossing the boundary of the continuum.Figure 4(b) presents the same plot but for slightly larger in-plane momentum q∥ =(0.05π, 0). The ω+ mode exhibits a qz dependence similar to that in Fig. 4(a). The ω−mode, however, shows a clear cosinelike dispersion along the qz direction and is situatedentirely above the continuum.In both Figs. 4(a) and (b), the intensity of the ω± modes displays a characteristic qzdependence, with “nodes” where the spectral weight almost vanishes at specific qz values.For example, the ω+ mode loses intensity around qz = 7π. These positions are parameterdependent, especially on the ratio of d/c, though the vanishing intensity of the ω− mode atqz = 0 is a robust feature. It is also important to note that the intensity of both ω± modesdoes not have 2π periodicity along the qz direction.For completeness, we also investigate the case of a large intrabilayer hopping tz = 0.3t.The spectral maps for various qz are shown in Fig. 5. Similar to the tz = 0.1t case, only asingle mode is present at qz = 0 [Fig. 5(a)], and its vanishing spectral weight at q∥ = (0, 0)identifies it as the ω+ mode. For qz ̸= 0 [Figs. 5(b)–(d)], two modes appear. A remarkablefeature in this large tz case is that the upper mode forms an upward-convex shape aroundq∥ = (0, 0) whereas the other mode is gapless and sharply defined even inside the continuum.At qz = 2π [Fig. 5(e)], two gapped modes are present. However, a key difference from thetz = 0.1t case is that the lower-energy mode in Fig. 5(e) is the one with zero spectral weightat q∥ = (0, 0), identifying it as the ω+ mode. Given the continuity of the modes with varyingqz, this suggests a reversal in the energy hierarchy of the ω± modes for a large tz. Specificallythe ω+ mode becomes gapless for qz ̸= 2nπ, while ω− mode is always gapped except for itsvanishing at qz = 0. This inversion can be traced back to the competing energy scales setby the LRC—through Vc and α—and by the intrabilayer hopping tz. When tz becomessufficiently large, the energy scale of the ω− mode can exceed the optical plasmon energy atqz = 2nπ, producing the observed reversal between the ω+ and ω− modes. This situationcan be realized for a system with a smaller intrabilayer distance. Designing materials in thisway would be highly interesting.As we observed in Fig. 5, the ω− mode has higher energy than the ω+ mode for largetz. This is confirmed by the results in Figs. 6(a) and (b), where the higher-energy modeshows zero spectral weight at qz = 0, characteristic of the ω− mode. The qz dependence ofthe ω− mode for large tz is markedly different from the tz = 0.1t case, showing a dip rather1100.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(a) qz = 0 tz = 0.3t   tz’ = 000.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(b) qz = 0.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(c) qz = π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(d) qz = 1.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(e) qz = 2πFIG. 5. Intensity maps of log10 |Imχ(q, ω)| for a sequence of qz around a region of q∥ = (0, 0) fortz = 0.3t: (a) qz = 0, (b) qz = 0.5π, (c) qz = π, (d) qz = 1.5π, and (e) qz = 2π. The white dottedcurve denotes the upper boundary of the particle-hole continuum. It becomes zero at q∥ = (0, 0)and qz = 0 in (a).1200.20.40.60.81.00 1 2 3 4 5 6 7 8ω / tqz/π<-2>0log10|Imχ(q,ω)|(a) q|| = (0.02π, 0) tz = 0.3ttz’ = 000.20.40.60.81.00 1 2 3 4 5 6 7 8ω / tqz/π<-2>0log10|Imχ(q,ω)|(b) q|| = (0.05π, 0)FIG. 6. qz dependence of ω+ mode (lower energy) and ω− mode (higher energy) at (a) q∥ =(0.02π, 0) and (b) q∥ = (0.05π, 0) for tz = 0.3t. The corresponding results of tz = 0.1t are shownin Fig. 4. The white dotted line is the upper boundary of the continuum spectrum and exhibits asharp drop at qz = 0 because of the vanishing of the ω− mode there.than a peak at qz = 2nπ. In contrast, the low-energy ω+ mode exhibits a significant disper-sive feature, particularly at q∥ = (0.02π, 0), where its energy spans a wide range betweenω = 0.06–0.44t. This mode is sharply defined even when it lies within the continuum, aconsequence of the very low spectral weight of the continuum itself near q∥ = (0, 0). Atq∥ = (0.05π, 0) shown in Fig. 6(b), the ω+ mode is less dispersive, but shows a similar overallfeature.B. Charge excitation spectra with t′z ̸= 0The previous analysis has focused on the case of zero interbilayer hopping (t′z = 0),though the LRC is included. Here, we extend our study to investigate the effect of a finitet′z. We begin by assuming a moderate value of t′z = tz/2 = 0.05t and later present resultsfor a larger t′z to capture the behavior observed for a large tz in the previous subsection.We present results in a manner consistent with Figs. 3(a)–(e), as shown in Figs. 7(a)–(e).At qz = 0 [Fig. 7(a)], the spectrum is nearly identical to the t′z = 0 case in Fig. 3(a), withonly the ω+ mode present. This indicates that the effect of t′z is negligible at qz = 0. Asimilar minor effect is seen at qz = 2π [Fig. 7(e)], where the ω− mode exists and is shifted toa slightly higher energy due to the inclusion of t′z. The most significant change induced bya finite t′z is observed at qz ̸= 2nπ, where the ω− mode, which was gapless for t′z = 0, nowbecomes gapped at q∥ = (0, 0). Consequently, for any qz ̸= 0, we have two gapped plasmon1300.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.8>0log10|Imχ(q,ω)|(a) qz = 0 tz = 0.1t   tz’ = tz/200.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.8>0log10|Imχ(q,ω)|(b) qz = 0.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.8>0log10|Imχ(q,ω)|(c) qz = π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.8>0log10|Imχ(q,ω)|(d) qz = 1.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.8>0log10|Imχ(q,ω)|(e) qz = 2πFIG. 7. Intensity maps of log10 |Imχ(q, ω)| for a sequence of qz around a region of q∥ = (0, 0) fortz = 0.1t and t′z = tz/2: (a) qz = 0, (b) qz = 0.5π, (c) qz = π, (d) qz = 1.5π, and (e) qz = 2π.The white dotted curve denotes the upper boundary of the particle-hole continuum. There is nocontinuum spectrum at q∥ = (0, 0) and qz = 0.1400.20.40.60.80 1 2 3 4 5 6 7 8ω / tqz/π<−2>0log10|Imχ(q,ω)|(a) q|| = (0.02π, 0) tz = 0.1ttz’ = tz/200.20.40.60.80 1 2 3 4 5 6 7 8ω / tqz/π<−2>0log10|Imχ(q,ω)|(b) q|| = (0.05π, 0)FIG. 8. qz dependence of ω+ mode (higher energy) and ω− mode (lower energy) at (a) q∥ =(0.02π, 0) and (b) q∥ = (0.05π, 0) for tz = 0.1t and t′z = tz/2. The white dotted curve denotes theupper boundary of the particle-hole continuum. There is a sharp drop at qz = 0 because of thedisappearance of the ω− mode there.modes in the presence of t′z.The qz dependence of the ω± modes is presented in Figs. 8(a) and (b) for q∥ = (0.02π, 0)and (0.05π, 0), respectively. A comparison with Figs. 4(a) and (b) shows that the overallbehavior is preserved, with the primary difference being a shift of the ω− mode to higherenergy due to the finite t′z. While the spectral intensity does not exhibit 2π periodicity, thelocation of the “nodes”, where the intensity is strongly suppressed, remain largely unaffectedby the inclusion of t′z. This suggests that the effect of t′z is relatively weak on the intensitydependence of the plasmon modes.For completeness, we also study the case of a large hopping tz = 0.3t with t′z = 0.15t. Aswe have previously established, the energy hierarchy of the ω± modes is interchanged for alarge tz, and this behavior persists in the presence of t′z.q-ω maps for this case are shown in Figs. 9(a)–(e). At qz = 0 [Fig. 9(a)], only theω+ mode is present, with its spectral weight vanishing at q∥ = (0, 0) as a consequence ofcharge conservation. For qz ̸= 0, the particle-hole continuum gains spectral weight even atq∥ = (0, 0), and the ω+ mode is realized close to this upper boundary around q∥ = (0, 0).In contrast, the ω− mode has a higher energy and is located above the continuum. As seenin Fig. 9(b), the ω− mode has a relatively low-spectral weight near qz = 0 and forms anupward-convex shape centered at q∥ = (0, 0) and ω ≈ 0.65t. This mode gains more spectralweight as qz increases as shown in Figs. 9(c)–(e), and its dispersion shows a small dependenceon qz. The low-energy ω+ mode enters slightly the continuum around q = (0, 0) in Figs. 9(b)1500.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(a) qz = 0 tz = 0.3t   tz’ = tz/200.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(b) qz = 0.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(c) qz = π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(d) qz = 1.5π00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<-1.7>0log10|Imχ(q,ω)|(e) qz = 2πFIG. 9. Intensity maps of log10 |Imχ(q, ω)| in the presence of a large tz = 0.3t and t′z = tz/2. (a)–(e) q∥ dependence for a sequence of qz around a region of q∥ = (0, 0): (a) qz = 0, (b) qz = 0.5π,(c) qz = π, (d) qz = 1.5π, and (e) qz = 2π. The white dotted curve denotes the upper boundary ofthe particle-hole continuum—there is no continuum spectrum at q∥ = (0, 0) and qz = 0 in (a).1600.20.40.60.81.00 1 2 3 4 5 6 7 8ω / tqz/π<−2.4>0log10|Imχ(q,ω)|(a) q|| = (0.02π, 0) tz = 0.3ttz’ = tz/200.20.40.60.81.00 1 2 3 4 5 6 7 8ω / tqz/π<−2.4>0log10|Imχ(q,ω)|(b) q|| = (0.05π, 0)FIG. 10. qz dependence of ω+ mode (lower energy) and ω− mode (higher energy) at (a) q∥ =(0.02π, 0) and (b) q∥ = (0.05π, 0) for tz = 0.3t and t′z = tz/2. The white dotted curve denotes theupper boundary of the particle-hole continuum. There is a sharp drop at qz = 0 because of thedisappearance of the ω− mode there.and (d). Its presence, despite being within the continuum, is due to the low-spectral weightof the continuum around q∥ = (0, 0). At qz = 2π [Fig. 9(e)], the ω+ mode is pushed upslightly above the continuum and its spectral intensity at q∥ = (0, 0) vanishes.Finally, we examine the qz dependence of the spectral intensity at fixed in-plane momentain Figs. 10(a) and (b) for q∥ = (0.02π, 0) and (0.05π, 0), respectively. The ω+ mode iswell defined when it is located above the continuum, and is somewhat blurred within thecontinuum, but a sharp peak is discernible at qz = 2nπ, particularly in Fig. 10(b). The ω−mode, which vanishes at qz = 0, exhibits a less dispersive feature along the qz direction, witha small dip at qz = 2nπ (n ̸= 0). The strong suppression of the ω− mode around qz = 7π isa feature also observed in the previous subsection.C. Vc dependence of dispersive modesThe t-J-V model contains a short-range interaction of the J-term as well as from thelocal constraint, and thus we can study the interplay with the long-range interaction V .In the case of qz = 0, we have only the ω+ mode with a gap [Fig. 3(a)]. This modeeventually becomes a gapless mode in the limit of Vc → 0, forming the zero-sound mode,as shown in Fig. 11, although the mode along the qx direction is realized very close to theupper boundary of the continuum in Fig. 11(b).At qz = π, both ω+ and ω− modes are present, with the ω− mode being gapless as shown1700.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<−1.7>0log10|Imχ(q,ω)|(a) Vc = t tz = 0.1t   tz’ = 0qz = 000.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω / t(q, q)/π              (q, 0)/π<−1.7>0log10|Imχ(q,ω)|(b) Vc = 0 (zero−sound mode)FIG. 11. Forming the zero-sound mode in the limit of Vc → 0: (a) Vc = t and (b) Vc = 0. It is theω+ mode that changes into the zero-sound mode.in Fig. 12(a). As the interaction strength Vc is reduced, the energy of the ω+ mode decreases,but it remains a gapped mode in the limit of Vc → 0. In contrast, the ω− mode retains itsgapless character, though the velocity is reduced as Vc decreases. While it might appearthat plasmon modes persist in the Vc = 0 limit, they are no longer plasmons. In particular,the ω− mode is realized inside the continuum as a peak of the continuum. On the otherhand, if we include t′z, both ω± modes are gapped even in the limit of Vc → 0.A unique feature is observed at qz = 2nπ. As previously shown in Fig. 3(e) both ω±modes are gapped when n ̸= 0. As Vc decreases [Fig. 12(b)], this gap is reduced. In contrastto the case of qz = 0 [Fig. 11(b)], the ω+ mode retains a gap in the limit of Vc → 0. Inthis limit, it is the ω− mode that changes into a gapless mode, namely evolves into thezero-sound-wavelike mode. This feature is the same even if there is a finite t′z.Corresponding spectra showing the qz dependence at q∥ = (0.05π, 0) are presented inFig. 13 in the limit of Vc = 0, where there are two branches. The analysis of this figureprovides three important insights. First, while the lower-energy branch at qz = 0 is the zero-sound mode, it smoothly connects to a branch at a finite qz. Second, because the V -termis the primary source of three-dimensional coupling, the system effectively reduces to beingtwo-dimensional in the Vc = 0 limit, resulting in a very weak dependence on qz. Third,1800.20.40.60.80 0.1 0.2 0.3ωp / t(q, 0)/πVc = 40tVc = 20tVc = 0(a) qz = π   tz = 0.1t   tz’ = 0   Γ = 10−3t00.20.40.60.80 0.1 0.2 0.3ωp / t(q, 0)/π(b) qz = 2πFIG. 12. Dispersive modes along the (0,0)-(0.3π, 0) direction for a sequence of values of Vc. Wehere employed a smaller broadening Γ = 10−3t to trace the dispersion precisely.00.10.20.30.40 1 2 3 4 5 6 7 8ω / tqz/π<−1.5>0log10|Imχ(q,ω)|Vc = 0tz = 0.1t   tz’ = 0   q|| = (0.05π, 0)FIG. 13. qz dependence of spectral weight for Vc = 0. While the zero-sound mode is well definedat qz = 0, it evolves smoothly into a lower-energy branch at qz ̸= 0. The higher-energy mode isabsent at qz = 0, but is realized at a finite qz.interestingly, the spectral intensity still exhibits a strong qz dependence due to the kinetichopping term along the z direction. Therefore, depending on the value of qz, the modes maybe detected as only one branch, although in principle two branches exist.19IV. DISCUSSIONSDespite the fundamental difference between our large-N theoretical framework and theRPA [39], we have found that our results for charge dynamics are strikingly similar in manyrespects, with a few notable differences.We have used a parameter set closely matching that of Ref. [39], with the exception of theCoulomb potential Vc and the anisotropic parameter α, and temperature. A key differencelies in our use of a significantly smaller value of Vc (in units of t). This disparity reflects thestrong correlation effects inherent in the t-J model, which lead to a notable band narrowingeffect. This strong correlation is a predominant factor, enabling the plasmon modes to existacross a wide range of momentum space as shown in Fig. 2(a). This sharply contrasts withthe RPA results [Fig. 2(b)], where plasmon modes are typically confined to the vicinity ofq∥ = (0, 0) and becomes heavily damped into the particle-hole continuum as they movefurther away.So far there is only one experimental RIXS report [38] in bilayer cuprates where a singleqy-scan was presented. Similar to weak-coupling calculations [39, 40], the present theory—amore adequate theory for cuprates—also suggests that the observed mode is the ω− modeas shown in Fig. 14. Because of a finite t′z, the ω− mode will be gapped at q∥ = (0, 0). To00.20.40.60.81.00.3 0.2 0.1 0 0.1 0.2 0.3ω [eV](q, q)/π              (0, q)/π<-1.7>0log10|Imχ(q,ω)|qz = 1.8πVc = 60eV   α = 13tz = 0.2t   tz’ = 0.05tFIG. 14. Comparison with the plasmon energy (solid squares) reported in Y-based cuprate su-perconductors in Ref. [38]. The experimental data are superimposed on the intensity map oflog10 |Imχ(q, ω)| computed at qz = 1.8π by using tz = 0.2t, t′z = 0.05t, Vc = 60 eV, and α = 13(corresponding to ϵ∥ = 1.68ϵ0 and ϵ⊥ = 1.16ϵ0), keeping the other parameters unchanged; t/2 isassumed to be 400 meV. The white dotted curve denotes the upper boundary of the particle-holecontinuum.20reinforce this conclusion, we call for more comprehensive data such as qx-, qy-, and qz-scans.In the bilayer t-J model two collective modes (zero-sound-wavelike modes) remain dueto electron correlations. In Appendix B we investigate whether these modes alone couldaccount for the experimental dispersion. Our analysis indicates that this scenario is lessconsistent with the data, especially regarding the slope of the measured dispersion, and thatthe plasmon interpretation offers a more coherent description of the observed features.V. CONCLUSIONSIn this work, we have constructed a strong-coupling theory of bilayer plasmons by em-ploying a large-N formalism for the t-J-V model. Our computational approach of chargeexcitation spectra was conducted in a matter that allows for a direct comparison with arecent RPA study [39]. Despite the fundamental differences in the theoretical framework,we have found a striking similarity in the plasmon dispersion and intensity maps. Thisagreement, however, is not without crucial distinctions. In our theory, the strong correla-tion effects inherent in the t-J model lead to a significant band narrowing, which in turnallows plasmon modes to remain well-defined across the entire Brillouin zone. This standsin sharp contrast to weak-coupling RPA calculations, where plasmons are typically heavilydamped away from the zone center.A unique contribution of our model is the insight gained by systematically varying theCoulomb interaction Vc. We have shown that as Vc is decreased, the plasmon modes changeinto two distinct modes. When qz = 2nπ (n ̸= 0), one of these modes remains gappedwhile the other becomes a gapless mode. When qz = 0, only the gapless mode—zero-soundmode—is present. We have also found that while the spectral intensity of these modes showsa strong dependence on qz, the mode energy itself is remarkably independent of qz.The ability of our strong-coupling theory to reproduce plasmon modes provides a compre-hensive framework for interpreting experimental data. Our results may offer an explanationfor the RIXS data recently obtained in Y-based cuprate superconductors. The ultimate testof the present strong-coupling theory will be the acquisition of more comprehensive RIXSdata in bilayer cuprates. If our model provides a coherent explanation for these future re-sults, it would offer compelling evidence that a strong-coupling approach is necessary fordescribing charge dynamics in these correlated systems.21ACKNOWLEDGMENTSThe authors thank M. Hepting and B. Keimer for valuable discussions about their data,and L. Benfatto, F. Gonzalez, W. Metzner, and I. Pomponio for illuminating discussions.H.Y. was supported by World Premier International Research Center Initiative (WPI),MEXT, Japan. M.B. and A.G. are indebted to warm hospitality of Max-Planck-Institute forSolid State Research. M.B. also thanks MANA Short-Term Invitation Program and warmhospitality in NIMS. A part of the results presented in this work was obtained by usingthe facilities of the CCT-Rosario Computational Center, member of the High PerformanceComputing National System (SNCAD, MincyT- Argentina).Appendix A: Large-N formalism of the bilayer t-J-V modelA major challenge in handling the t-J model arises from the non-double-occupancy con-straint. To rigorously enforce the local constraint [Eq. (2)], we express the Hamiltonian interms of the Hubbard operators X̂[57]. The constraint is then implicitly described by the al-gebra of these operators: c̃†iσ,α = X̂σ0iα , c̃iσ,α = X̂0σiα , S+i,α = X̂↑↓iα , S−i,α = X̂↓↑iα , ni,α = X̂↑↑iα +X̂↓↓iα ,and X̂00iα describes the number of doped holes. The z component of the spin operator isSzi,α = 12(X̂↑↑iα − X̂↓↓iα ). The operators X̂σ0iα and X̂0σiα are called fermionlike, whereas theoperators X̂σσ′iα and X̂00iα are bosonlike.The Hamiltonian in Eq. (1) can be expressed in terms of the Hubbard operators as:H =−∑i,j,σ,α,βtαβij Xσ0iαX0σjβ − µ∑i,σ,αXσσiα +J2∑⟨i,j⟩,σ,σ′,α(Xσσ′iα Xσ′σjα −Xσσiα Xσ′σ′jα)+J⊥4∑i,σ,σ′,α̸=β(Xσσ′iα Xσ′σiβ −Xσσiα Xσ′σ′iβ)+12∑i̸=j,σ,σ′,α,βV αβij Xσσiα Xσ′σ′jβ . (A1)The formalism starts with the construction of a first-order classical Lagrangian using theFaddeev-Jackiw and Dirac methods [58–60]. In this representation, the fermionlike (boson-like) Hubbard operators are associated with Grassmann (usual bosonic) variables. Next, alarge-N expansion is applied to the spin projection, extending it from σ =↑, ↓ to p = 1, ..., Nand rescaling the amplitude as tαβij /N , J/N , J⊥/N , and V αβij /N to ensure a finite theory inthe limit of N → ∞. Using the conditionXpp′iα =Xp0iαX0p′iαX00iα(A2)22for J- and J⊥-terms, we write the fermionlike fields asf †ip,α =1√Nδ/2Xp0iα , (A3)fip,α =1√Nδ/2X0piα , (A4)where δ is the hole doping away from half-filling. Xσσiα in J-, J⊥- and V αβij -terms are, on theother hand, treated by utilizing the local constraint X00iα +∑p Xppiα = N/2 and this constraintis imposed by introducing the Lagrange multiplier λiα.The fields X00iα and λiα are expressed in terms of their static mean-field components anddynamic fluctuations:X00iα = Nδ2(1 + δRiα) , (A5)λiα = λ0 + δλiα , (A6)where δRiα denotes the fluctuation of the hole density at site i in layer α; δλiα is thefluctuation of the Lagrange multiplier to enforce the constraint against double occupancy.The resulting effective Lagrangian includes two distinct four-fermion interaction terms,one from the in-plane exchange interaction J , and the other from the out-of-plane interactionJ⊥. To decouple these terms, we introduce Hubbard-Stratonovich fields ∆ij,α and ∆′i:∆ij,α =J2∑pf †jp,αfip,α√(1 + δRiα)(1 + δRjα), (A7)∆′i =J⊥4∑pf †ip,1fip,2√(1 + δRi1)(1 + δRi2). (A8)The fields ∆ij,α and ∆′i describe bond-charge fluctuations in the intralayer and intrabilayer,respectively. Since i and j are nearest-neighbor sites on the square lattice, we may write∆ij,α = ∆ηiα where η = x or y. We parametrize those fields as:∆ηiα = χ (1 + rηiα + iAηiα) , (A9)∆′i = χ′ (1 + r⊥,i + iA⊥,i) , (A10)where rηiα and Aηiα (r⊥,i and A⊥,i) represent the real and imaginary parts of the in-plane(out-of-plane) bond-field fluctuations, respectively, and χ (χ′) is the corresponding staticmean-field value.23Finally, the terms involving 1/√1 + δRiα are expanded perturbatively in powers of δRiα.This expansion systematically organizes the interactions in powers of 1/N , thus controllingthe hierarchy of contributions in the large-N formalism. The effective theory of Eq. (1) isthen described in terms of fermions, bosons, and their interactions.In the large-N formalism for bilayer systems, we introduce a 14-component bosonic fieldas a basis:δXa = (δR1, δλ1, δR2, δλ2, rx1 , ry1 , Ax1 , Ay1, rx2 , ry2 , Ax2 , Ay2, r⊥, A⊥) , (A11)where the site index is omitted for clarity.Following Refs. [52, 61, 62], the Feynman rules applied to the effective theory give a14× 14 bare bosonic propagator D(0)ab (q, iωn) that is O(1/N):[D(0)ab (q, iωn)]−1= ND(0)A D(0)B 0 0 0D(0)∗B D(0)A 0 0 00 0 D(0)C 0 00 0 0 D(0)C 00 0 0 0 D(0)D, (A12)where ωn is a bosonic Matsubara frequency and the matrices D(0)A-D are:D(0)A = δ22[V (q)2− J(q)]δ2δ20 , (A13)D(0)B = δ22[V ′(q)2− J ′(q)]00 0 , (A14)D(0)C =4χ2J0 0 00 4χ2J0 00 0 4χ2J00 0 0 4χ2J , (A15)D(0)D = 4χ′2J⊥00 4χ′2J⊥ . (A16)24k,�k',�q,a���,a =k,�k',�q,a���,ab =q',bb q ak� �G(0) =�� D(0) = ab ab =b a+b aD-1 =                   = [D(0)]-1 - ab abab (             )-1(a)(b)(c)(d)FIG. 15. (a) Bare fermionic G(0)αβ and bosonic D(0)ab propagators, solid and dashed lines, respectively.(b) Λαβ,a and Λαβ,ab are the three- and four-legs vertices, respectively. (c) Πab is the bosonic self-energy. (d) Double dashed line is the dressed bosonic propagator Dab.The large-N formalism yields the bare bosonic propagator as a 14 × 14 matrix [seeEq. (A12)]. However, because we focus on on-site charge excitations, including plasmons[63], we restrict our analysis to the corresponding 4 × 4 submatrix of the bare bosonicpropagator given by[D(0)ab (q, iωn)]−1= Nδ22[V (q)2− J(q)]δ2δ22[V ′(q)2− J ′(q)]0δ20 0 0δ22[V ′∗(q)2− J ′∗(q)]0 δ22[V (q)2− J(q)]δ20 0 δ20 , (A17)where the matrix elements are limited to a, b = δRa, δλa.Interactions between bosons and fermions are governed by three- and four-leg vertices.The three-leg vertex Λαβ,a represents the interaction of a fermion from layer α that ends atlayer β after interacting with a boson δXa. We can write this vertex asΛαβ,a(k, k′, q) = Λ̃αβ,a(k, k′, q)ei(−kzdα+k′zdβ+qzda)δ(k − k′ − q) (A18)where k ≡ k, iνn, q ≡ q, iωn and da is equal to d1 or d2 depending on the plane of the bosonδXa, with d1 = 0 and d2 = d/c. After executing δ(k − k′ − q), the nonzero components ofΛ̃αβ,a(k, k′, q) in the charge sector become25Λ̃αα,a(k, q) = −[iνn −iωn2+ µ+ 2χ∑η=x,ycos(kη −qη2)cosqη2, 1](A19)for each component a = δRα, δλα. For α ̸= βΛ̃αβ,a(k, q) =−(χ′2,χ′2)(A20)for each component a = δR1, δR2.The four-leg vertex Λαβ,ab represents a fermion from layer α that ends in layer β afterinteracting with the bosons δXa and δXb. We write this vertex asΛαβ,ab(k, k′, q, q′) = Λ̃αβ,ab(k, k′, q, q′)ei(−kzdα+k′zdβ+qzda+q′zdb)δ(k − k′ − q − q′) , (A21)where, after executing δ(k− k′ − q− q′), the nonzero components of Λ̃αβ,ab(k, k′, q, q′) in thecharge sector becomeΛ̃αα,ab(k, q, q′) = Fk,q,q′ 1/21/2 0 (A22)for each component a, b = δRα, δλα, whereFk,q,q′ = iνn −iωn + iω′n2+ µ+ χ∑η=x,ycos(kη −qη + q′η2)×[cos(qη + q′η2)+ cosqη2cosq′η2]. (A23)For α ̸= βΛ̃αβ,ab =χ′8 3 11 3 (A24)for each component a, b = δR1, δR2. Note that the vertices are O(1).Using the propagators [Eq. (3) and Fig. 15(a)] and vertices [Fig. 15(b)] the bosonic self-energy Πab(q, iωn) [Fig. 15(c)] is computed, considering both Hartree and bubble diagrams.The Dyson equation [Fig. 15(d)] yields the dressed bosonic propagator Dab(q, iωn) asD−1ab (q, iωn) =[D(0)ab (q, iωn)]−1− Πab(q, iωn) . (A25)26Focusing on the 4× 4 charge sector of the self-energy, the analytical expressions of eachcomponent are:Π11(q, iωn) =− N16Ns2∑α,β=1∑k[nF(εαk−q)− nF(εβk)](ε̃βk − ε̃αk−q)+(ε̃βk + ε̃αk−q)2gαβ ,(A26)Π12(q, iωn) =− N8Ns2∑α,β=1∑k(ε̃βk + ε̃αk−q)gαβ , (A27)Π13(q, iωn) =− e−iqzdcN16Ns2∑α,β=1(−1)α+β∑kε⊥∗k−qε⊥k|ε⊥k−q||ε⊥k |{[nF(εαk−q)− nF(εβk)](ε̃βk − ε̃αk−q)+(ε̃βk + ε̃αk−q)2gαβ}, (A28)Π31(q, iωn) =− eiqzdcN16Ns2∑α,β=1(−1)α+β∑kε⊥k−qε⊥∗k|ε⊥k−q||ε⊥k |{[nF(εαk−q)− nF(εβk)](ε̃βk − ε̃αk−q)+(ε̃βk + ε̃αk−q)2gαβ}, (A29)Π14(q, iωn) =− e−iqzdcN8Ns2∑α,β=1(−1)α+β∑kε⊥∗k−qε⊥k|ε⊥k−q||ε⊥k |(ε̃βk + ε̃αk−q)gαβ (A30)Π41(q, iωn) =− eiqzdcN8Ns2∑α,β=1(−1)α+β∑kε⊥k−qε⊥∗k|ε⊥k−q||ε⊥k |(ε̃βk + ε̃αk−q)gαβ (A31)Π22(q, iωn) =− N4Ns2∑α,β=1∑kgαβ (A32)Π24(q, iωn) =− e−iqzdcN4Ns2∑α,β=1(−1)α+β∑kε⊥∗k−qε⊥k|ε⊥k−q||ε⊥k |gαβ , (A33)Π42(q, iωn) =− eiqzdcN4Ns2∑α,β=1(−1)α+β∑kε⊥k−qε⊥∗k|ε⊥k−q||ε⊥k |gαβ , (A34)Π21(q, iωn) =Π34(q, iωn) = Π43(q, iωn) = Π12(q, iωn) , (A35)Π23(q, iωn) =Π14(q, iωn) , (A36)Π32(q, iωn) =Π41(q, iωn) , (A37)Π33(q, iωn) =Π11(q, iωn) , (A38)Π44(q, iωn) =Π22(q, iωn) , (A39)27where ε̃αk is equal to εαk with µ = χ = χ′ = 0, andgαβ =nF(εαk−q)− nF(εβk)iωn + εαk−q − εβk. (A40)Appendix B: Analysis of experimental data with Vc = 0In Sec. III C we showed that two collective modes persist even when the long-rangeCoulomb interaction is switched off (Vc = 0). In this section we examine whether theexperimental dispersion in Fig. 14 could be interpreted within this scenario.A first qualitative consideration already points to a limitation: the lower-energy zero-sound-wavelike mode necessarily decreases linearly to zero as q∥ → (0, 0), whereas theexperimental data show no indication of such behavior. This implies that, if a zero-sound-based interpretation were to be considered, the relevant branch would need to be the higher-energy mode.Figure 16(a) shows the intensity map for Vc = 0 using the same band parameters thatsuccessfully describe the plasmon dispersion in Fig. 14. The resulting zero-sound-wavelikebranches lie at substantially lower energies than the experimental data, indicating thatthese parameters cannot account for the measured spectrum. To counteract this effect, we00.20.40.60.81.00 0.05 0.10 0.15 0.20 0.25 0.30ω [eV](0, q)/π<-1>0log10|Imχ(q,ω)|qz = 1.8πVc = 0(a) tz = 0.2t   tz’ = 0.05t00.20.40.60.81.00 0.05 0.10 0.15 0.20 0.25 0.30ω [eV](0, q)/π<-1>0log10|Imχ(q,ω)|(b) tz = 0.4t   tz’ = 0.15tFIG. 16. Comparison with the plasmon energy (solid squares) reported in Y-based cuprate su-perconductors in Ref. [38]. The experimental data are superimposed on the intensity map oflog10 |Imχ(q, ω)| computed at qz = 1.8π with Vc = 0 by using (a) tz = 0.2t, t′z = 0.05t (same asin Fig. 14), and (b) tz = 0.4t, t′z = 0.15t, while keeping the other parameters unchanged; t/2 isassumed to be 400 meV. The white dotted curve denotes the upper boundary of the particle-holecontinuum.28consider larger values of tz and t′z and display the corresponding results in Fig. 16(b). Whilethe overall energy scale then becomes comparable to the measured one, the slope of thecomputed dispersion remains inconsistent with the observed trend. Moreover, the requiredvalues of tz and t′z appear too large for bilayer cuprates, suggesting that a description basedsolely on zero-sound-wavelike modes is unlikely.These analyses lead us to conclude that the experimental features cannot be satisfactorilydescribed by the collective modes of the pure bilayer t-J model. Instead, their dispersion andspectral characteristics are more naturally and consistently interpreted in terms of plasmons.[1] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, “From quantum matterto high-temperature superconductivity in copper oxides,” Nature 518, 179–186 (2015).[2] P. W. Anderson, “The resonating valence bond state in La2CuO4 and superconductivity,”Science 235, 1196–1198 (1987).[3] Luuk J. P. Ament, Michel van Veenendaal, Thomas P. Devereaux, John P. Hill, and Jeroenvan den Brink, “Resonant inelastic x-ray scattering studies of elementary excitations,” Rev.Mod. Phys. 83, 705–767 (2011).[4] Frank M. F. de Groot, Maurits W. Haverkort, Hebatalla Elnaggar, Amélie Juhin, Ke-JinZhou, and Pieter Glatzel, “Resonant inelastic x-ray scattering,” Nature Reviews MethodsPrimers 4, 45 (2024).[5] G. Ghiringhelli, M. Le Tacon, M. Minola, S. Blanco-Canosa, C. Mazzoli, N. B. Brookes, G. M.De Luca, A. Frano, D. G. Hawthorn, F. He, T. Loew, M. Moretti Sala, D. C. Peets, M. Salluzzo,E. Schierle, R. Sutarto, G. A. Sawatzky, E. Weschke, B. Keimer, and L. Braicovich, “Long-Range Incommensurate Charge Fluctuations in (Y,Nd)Ba2Cu3O6+x,” Science 337, 821–825(2012).[6] J. Chang, E. Blackburn, A. T. Holmes, N. B. Christensen, J. Larsen, J. Mesot, Ruixing Liang,D. A. Bonn, W. N. Hardy, A. Watenphul, M. v. Zimmermann, E. M. Forgan, and S. M.Hayden, “Direct observation of competition between superconductivity and charge densitywave order in YBa2Cu3O6.67,” Nat. Phys. 8, 871 (2012).[7] A. J. Achkar, R. Sutarto, X. Mao, F. He, A. Frano, S. Blanco-Canosa, M. Le Tacon, G. Ghir-inghelli, L. Braicovich, M. Minola, M. Moretti Sala, C. Mazzoli, Ruixing Liang, D. A. Bonn,29http://dx.doi.org/ 10.1038/nature14165http://dx.doi.org/10.1126/science.235.4793.1196http://dx.doi.org/ 10.1038/s43586-024-00322-6http://dx.doi.org/ 10.1038/s43586-024-00322-6http://dx.doi.org/10.1126/science.1223532http://dx.doi.org/10.1126/science.1223532W. N. Hardy, B. Keimer, G. A. Sawatzky, and D. G. Hawthorn, “Distinct Charge Ordersin the Planes and Chains of Ortho-III-Ordered YBa2Cu3O6+δ Superconductors Identified byResonant Elastic X-ray Scattering,” Phys. Rev. Lett. 109, 167001 (2012).[8] E. Blackburn, J. Chang, M. Hücker, A. T. Holmes, N. B. Christensen, Ruixing Liang, D. A.Bonn, W. N. Hardy, U. Rütt, O. Gutowski, M. v. Zimmermann, E. M. Forgan, and S. M.Hayden, “X-Ray Diffraction Observations of a Charge-Density-Wave Order in Superconduct-ing Ortho-II YBa2Cu3O6.54 Single Crystals in Zero Magnetic Field,” Phys. Rev. Lett. 110,137004 (2013).[9] S. Blanco-Canosa, A. Frano, E. Schierle, J. Porras, T. Loew, M. Minola, M. Bluschke,E. Weschke, B. Keimer, and M. Le Tacon, “Resonant x-ray scattering study of charge-densitywave correlations in YBa2Cu3O6+x,” Phys. Rev. B 90, 054513 (2014).[10] R. Comin, A. Frano, M. M. Yee, Y. Yoshida, H. Eisaki, E. Schierle, E. Weschke, R. Sutarto,F. He, A. Soumyanarayanan, Yang He, M. Le Tacon, I. S. Elfimov, Jennifer E. Hoffman, G. A.Sawatzky, B. Keimer, and A. Damascelli, “Charge Order Driven by Fermi-Arc Instability inBi2Sr2−xLaxCuO6+δ,” Science 343, 390–392 (2014).[11] W. Tabis, Y. Li, M. Le Tacon, L. Braicovich, A. Kreyssig, M. Minola, G. Dellea, E. Weschke,M. J. Veit, M. Ramazanoglu, A. I. Goldman, T. Schmitt, G. Ghiringhelli, N. Barisic, M. K.Chan, C. J. Dorow, G. Yu, X. Zhao, B. Keimer, and M. Greven, “Charge order and itsconnection with Fermi-liquid charge transport in a pristine high-Tc cuprate,” Nat. Commun.5, 5875 (2014).[12] Eduardo H. da Silva Neto, Pegor Aynajian, Alex Frano, Riccardo Comin, Enrico Schierle,Eugen Weschke, András Gyenis, Jinsheng Wen, John Schneeloch, Zhijun Xu, Shimpei Ono,Genda Gu, Mathieu Le Tacon, and Ali Yazdani, “Ubiquitous interplay between charge or-dering and high-temperature superconductivity in cuprates,” Science 343, 393–396 (2014).[13] M. Hashimoto, G. Ghiringhelli, W.-S. Lee, G. Dellea, A. Amorese, C. Mazzoli, K. Kum-mer, N. B. Brookes, B. Moritz, Y. Yoshida, H. Eisaki, Z. Hussain, T. P. Devereaux, Z.-X.Shen, and L. Braicovich, “Direct observation of bulk charge modulations in optimally dopedBi1.5Pb0.6Sr1.54CaCu2O8+δ,” Phys. Rev. B 89, 220511 (2014).[14] Y. Y. Peng, M. Salluzzo, X. Sun, A. Ponti, D. Betto, A. M. Ferretti, F. Fumagalli, K. Kummer,M. Le Tacon, X. J. Zhou, N. B. Brookes, L. Braicovich, and G. Ghiringhelli, “Direct obser-vation of charge order in underdoped and optimally doped Bi2(Sr,La)2CuO6+δ by resonant30http://dx.doi.org/ 10.1103/PhysRevB.90.054513http://dx.doi.org/ 10.1126/science.1242996http://dx.doi.org/10.1038/ncomms6875http://dx.doi.org/10.1038/ncomms6875http://dx.doi.org/ 10.1126/science.1243479http://dx.doi.org/10.1103/PhysRevB.89.220511inelastic x-ray scattering,” Phys. Rev. B 94, 184511 (2016).[15] L. Chaix, G. Ghiringhelli, Y. Y. Peng, M. Hashimoto, B. Moritz, K. Kummer, N. B. Brookes,Y. He, S. Chen, S. Ishida, Y. Yoshida, H. Eisaki, M. Salluzzo, L. Braicovich, Z. X. Shen, T. P.Devereaux, and W. S. Lee, “Dispersive charge density wave excitations in Bi2Sr2CaCu2O8+δ,”Nat. Phys. 13, 952 (2017).[16] R. Arpaia, S. Caprara, R. Fumagalli, G. De Vecchi, Y. Y. Peng, E. Andersson, D. Betto,G. M. De Luca, N. B. Brookes, F. Lombardi, M. Salluzzo, L. Braicovich, C. Di Castro,M. Grilli, and G. Ghiringhelli, “Dynamical charge density fluctuations pervading the phasediagram of a Cu-based high-Tc superconductor,” Science 365, 906–910 (2019).[17] B. Yu, W. Tabis, I. Bialo, F. Yakhou, N. B. Brookes, Z. Anderson, Y. Tang, G. Yu, andM. Greven, “Unusual dynamic charge correlations in simple-tetragonal HgBa2CuO4+δ,” Phys.Rev. X 10, 021059 (2020).[18] W. S. Lee, Ke-Jin Zhou, M. Hepting, J. Li, A. Nag, A. C. Walters, M. Garcia-Fernandez, H. C.Robarts, M. Hashimoto, H. Lu, B. Nosarzewski, D. Song, H. Eisaki, Z. X. Shen, B. Moritz,J. Zaanen, and T. P. Devereaux, “Spectroscopic fingerprint of charge order melting driven byquantum fluctuations in a cuprate,” Nat. Phys. 17, 53–57 (2021).[19] Haiyu Lu, Makoto Hashimoto, Su-Di Chen, Shigeyuki Ishida, Dongjoon Song, Hiroshi Eisaki,Abhishek Nag, Mirian Garcia-Fernandez, Riccardo Arpaia, Giacomo Ghiringhelli, LucioBraicovich, Jan Zaanen, Brian Moritz, Kurt Kummer, Nicholas B. Brookes, Ke-Jin Zhou,Zhi-Xun Shen, Thomas P. Devereaux, and Wei-Sheng Lee, “Identification of a characteristicdoping for charge order phenomena in Bi-2212 cuprates via RIXS,” Phys. Rev. B 106, 155109(2022).[20] Riccardo Arpaia, Leonardo Martinelli, Marco Moretti Sala, Sergio Caprara, Abhishek Nag,Nicholas B. Brookes, Pietro Camisa, Qizhi Li, Qiang Gao, Xingjiang Zhou, Mirian Garcia-Fernandez, Ke-Jin Zhou, Enrico Schierle, Thilo Bauch, Ying Ying Peng, Carlo Di Castro,Marco Grilli, Floriana Lombardi, Lucio Braicovich, and Giacomo Ghiringhelli, “Signatureof quantum criticality in cuprates by charge density fluctuations,” Nat. Commun. 14, 7198(2023).[21] (), the charge order observed in La-based cuprates is often discussed in terms of the spin-charge stripe order [64]. It is an interesting subject to investigate its possible connection tothe newly observed charge ordering tendency around qx ∼ 0.6π in other hole-doped cuprates.31http://dx.doi.org/ 10.1103/PhysRevB.94.184511[22] Eduardo H. da Silva Neto, Riccardo Comin, Feizhou He, Ronny Sutarto, Yeping Jiang,Richard L. Greene, George A. Sawatzky, and Andrea Damascelli, “Charge ordering in theelectron-doped superconductor Nd2−xCexCuO4,” Science 347, 282–285 (2015).[23] Eduardo H. da Silva Neto, Biqiong Yu, Matteo Minola, Ronny Sutarto, Enrico Schierle, FabioBoschini, Marta Zonno, Martin Bluschke, Joshua Higgins, Yangmu Li, Guichuan Yu, Eu-gen Weschke, Feizhou He, Mathieu Le Tacon, Richard L. Greene, Martin Greven, George A.Sawatzky, Bernhard Keimer, and Andrea Damascelli, “Doping-dependent charge order cor-relations in electron-doped cuprates,” Science Advances 2, e1600782 (2016).[24] E. H. da Silva Neto, M. Minola, B. Yu, W. Tabis, M. Bluschke, D. Unruh, H. Suzuki, Y. Li,G. Yu, D. Betto, K. Kummer, F. Yakhou, N. B. Brookes, M. Le Tacon, M. Greven, B. Keimer,and A. Damascelli, “Coupling between dynamic magnetic and charge-order correlations in thecuprate superconductor Nd2−xCexCuO4,” Phys. Rev. B 98, 161114 (2018).[25] K. Ishii, K. Tsutsui, Y. Endoh, T. Tohyama, S. Maekawa, M. Hoesch, K. Kuzushita, M. Tsub-ota, T. Inami, J. Mizuki, Y. Murakami, and K. Yamada, “Momentum Dependence of ChargeExcitations in the Electron-Doped Superconductor Nd1.85Ce0.15CuO4: A Resonant InelasticX-Ray Scattering Study,” Phys. Rev. Lett. 94, 207003 (2005).[26] K. Ishii, M. Fujita, T. Sasaki, M. Minola, G. Dellea, C. Mazzoli, K. Kummer, G. Ghiringhelli,L. Braicovich, T. Tohyama, K. Tsutsumi, K. Sato, R. Kajimoto, K. Ikeuchi, K. Yamada,M. Yoshida, M. Kurooka, and J. Mizuki, “High-energy spin and charge excitations in electron-doped copper oxide superconductors,” Nat. Commun. 5, 3714 (2014).[27] W. S. Lee, J. J. Lee, E. A. Nowadnick, S. Gerber, W. Tabis, S. W. Huang, V. N. Strocov,E. M. Motoyama, G. Yu, B. Moritz, H. Y. Huang, R. P. Wang, Y. B. Huang, W. B. Wu, C. T.Chen, D. J. Huang, M. Greven, T. Schmitt, Z. X. Shen, and T. P. Devereaux, “Asymmetryof collective excitations in electron- and hole-doped cuprate superconductors,” Nat. Phys. 10,883–889 (2014).[28] Andrés Greco, Hiroyuki Yamase, and Mat́ıas Bejas, “Plasmon excitations in layered high-Tccuprates,” Phys. Rev. B 94, 075139 (2016).[29] Kenji Ishii, Takami Tohyama, Shun Asano, Kentaro Sato, Masaki Fujita, Shuichi Wakimoto,Kenji Tustsui, Shigetoshi Sota, Jun Miyawaki, Hideharu Niwa, Yoshihisa Harada, JonathanPelliciari, Yaobo Huang, Thorsten Schmitt, Yoshiya Yamamoto, and Jun’ichiro Mizuki, “Ob-servation of momentum-dependent charge excitations in hole-doped cuprates using resonant32http://dx.doi.org/ 10.1126/science.1256441http://dx.doi.org/ 10.1126/sciadv.1600782http://dx.doi.org/10.1103/PhysRevB.98.161114http://dx.doi.org/ 10.1103/PhysRevLett.94.207003http://dx.doi.org/ http://dx.doi.org/10.1038/ncomms4714 10.1038/ncomms4714http://dx.doi.org/http://dx.doi.org/10.1038/nphys3117 10.1038/nphys3117http://dx.doi.org/http://dx.doi.org/10.1038/nphys3117 10.1038/nphys3117http://dx.doi.org/10.1103/PhysRevB.94.075139inelastic x-ray scattering at the oxygen K edge,” Phys. Rev. B 96, 115148 (2017).[30] G. Dellea, M. Minola, A. Galdi, D. Di Castro, C. Aruta, N. B. Brookes, C. J. Jia, C. Mazzoli,M. Moretti Sala, B. Moritz, P. Orgiani, D. G. Schlom, A. Tebano, G. Balestrino, L. Braicovich,T. P. Devereaux, L. Maritato, and G. Ghiringhelli, “Spin and charge excitations in artificialhole- and electron-doped infinite layer cuprate superconductors,” Phys. Rev. B 96, 115117(2017).[31] M. Hepting, L. Chaix, E. W. Huang, R. Fumagalli, Y. Y. Peng, B. Moritz, K. Kummer, N. B.Brookes, W. C. Lee, M. Hashimoto, T. Sarkar, J.-F. He, C. R. Rotundu, Y. S. Lee, R. L.Greene, L. Braicovich, G. Ghiringhelli, Z. X. Shen, T. P. Devereaux, and W. S. Lee, “Three-dimensional collective charge excitations in electron-doped copper oxide superconductors,”Nature 563, 374–378 (2018).[32] Andrés Greco, Hiroyuki Yamase, and Mat́ıas Bejas, “Origin of high-energy charge excitationsobserved by resonant inelastic x-ray scattering in cuprate superconductors,” Commun. Phys.2, 3 (2019).[33] M. Hepting, M. Bejas, A. Nag, H. Yamase, N. Coppola, D. Betto, C. Falter, M. Garcia-Fernandez, S. Agrestini, Ke-Jin Zhou, M. Minola, C. Sacco, L. Maritato, P. Orgiani, H. I.Wei, K. M. Shen, D. G. Schlom, A. Galdi, A. Greco, and B. Keimer, “Gapped collectivecharge excitations and interlayer hopping in cuprate superconductors,” Phys. Rev. Lett. 129,047001 (2022).[34] Andrés Greco, Hiroyuki Yamase, and Mat́ıas Bejas, “Close inspection of plasmon excitationsin cuprate superconductors,” Phys. Rev. B 102, 024509 (2020).[35] Alexander L Fetter, “Electrodynamics of a layered electron gas. II. periodic array,” Annals ofPhysics 88, 1–25 (1974).[36] A. Griffin and A. J. Pindor, “Plasmon dispersion relations and the induced electron interactionin oxide superconductors: Numerical results,” Phys. Rev. B 39, 11503–11514 (1989).[37] Akira Iyo, Yasumoto Tanaka, Hijiri Kito, Yasuharu Kodama, Parasharam M. Shirage, DilipD. Shivagan, Hirofumi Matsuhata, Kazuyasu Tokiwa, and Tsuneo Watanabe, “Tc vs n Rela-tionship for Multilayered High-Tc Superconductors,” Journal of the Physical Society of Japan76, 094711 (2007).[38] M. Bejas, V. Zimmermann, D. Betto, T. D. Boyko, R. J. Green, T. Loew, N. B. Brookes,G. Cristiani, G. Logvenov, M. Minola, B. Keimer, H. Yamase, A. Greco, and M. Hepting,33http://dx.doi.org/10.1103/PhysRevB.96.115148http://dx.doi.org/10.1103/PhysRevB.96.115117http://dx.doi.org/10.1103/PhysRevB.96.115117http://dx.doi.org/10.1038/s41586-018-0648-3http://dx.doi.org/ 10.1038/s42005-018-0099-zhttp://dx.doi.org/ 10.1038/s42005-018-0099-zhttp://dx.doi.org/ 10.1103/PhysRevLett.129.047001http://dx.doi.org/ 10.1103/PhysRevLett.129.047001http://dx.doi.org/10.1103/PhysRevB.102.024509http://dx.doi.org/https://doi.org/10.1016/0003-4916(74)90397-2http://dx.doi.org/https://doi.org/10.1016/0003-4916(74)90397-2http://dx.doi.org/ 10.1103/PhysRevB.39.11503http://dx.doi.org/10.1143/JPSJ.76.094711http://dx.doi.org/10.1143/JPSJ.76.094711“Plasmon dispersion in bilayer cuprate superconductors,” Phys. Rev. B 109, 144516 (2024).[39] Hiroyuki Yamase, “Theory of charge dynamics in bilayer electron system with long-rangecoulomb interaction,” Phys. Rev. B 111, 085138 (2025).[40] Niccolò Sellati and Lara Benfatto, “Ghost josephson plasmon in bilayer superconductors,”Phys. Rev. B 111, 104509 (2025).[41] S. Nakata, M. Bejas, J. Okamoto, K. Yamamoto, D. Shiga, R. Takahashi, H. Y. Huang,H. Kumigashira, H. Wadati, J. Miyawaki, S. Ishida, H. Eisaki, A. Fujimori, A. Greco, H. Ya-mase, D. J. Huang, and H. Suzuki, “Out-of-phase plasmon excitations in the trilayer cuprateBi2Sr2Ca2Cu3O10+δ,” Phys. Rev. B 111, 165141 (2025).[42] M. Mitrano, A. A. Husain, S. Vig, A. Kogar, M. S. Rak, S. I. Rubeck, J. Schmalian, B. Uchoa,J. Schneeloch, R. Zhong, G. D. Gu, and P. Abbamonte, “Anomalous density fluctuations ina strange metal,” Proc. Natl. Acad. Sci. U. S. A. 115, 5392–5396 (2018).[43] Aurelio Romero-Bermúdez, Alexander Krikun, Koenraad Schalm, and Jan Zaanen, “Anoma-lous attenuation of plasmons in strange metals and holography,” Phys. Rev. B 99, 235149(2019).[44] S. T. Van den Eede, T. J. N. van Stralen, C. F. J. Flipse, and H. T. C. Stoof, “Plasmons ina layered strange metal using the gauge-gravity duality,” Phys. Rev. B 109, 085119 (2024).[45] H. Sun, M. Huo, X. Hu, J. Li, Z. Liu, Y. Han, L. Tang, Z. Mao, P. Yang, B. Wang, J. Cheng,D.-X. Yao, G.-M. Zhang, and M. Wang, “Signatures of superconductivity near 80 K in anickelate under high pressure,” Nature 621, 493–498 (2023).[46] J. Hou, P.-T. Yang, Z.-Y. Liu, J.-Y. Li, P.-F. Shan, L. Ma, G. Wang, N.-N. Wang, H.-Z. Guo,J.-P. Sun, Y. Uwatoko, M. Wang, G.-M. Zhang, B.-S. Wang, and J.-G. Cheng, “Emergence ofHigh-Temperature Superconducting Phase in Pressurized La3Ni2O7 Crystals,” Chinese Phys.Lett. 40, 117302 (2023).[47] Adriana Foussats and Andrés Greco, “Large-N expansion based on the Hubbard operatorpath integral representation and its application to the t-J model,” Phys. Rev. B 65, 195107(2002).[48] Abhishek Nag, M. Zhu, Mat́ıas Bejas, J. Li, H. C. Robarts, Hiroyuki Yamase, A. N. Petsch,D. Song, H. Eisaki, A. C. Walters, M. Garćıa-Fernández, Andrés Greco, S. M. Hayden, and Ke-Jin Zhou, “Detection of Acoustic Plasmons in Hole-Doped Lanthanum and Bismuth CuprateSuperconductors Using Resonant Inelastic X-Ray Scattering,” Phys. Rev. Lett. 125, 25700234http://dx.doi.org/10.1103/PhysRevB.109.144516http://dx.doi.org/10.1103/PhysRevB.111.085138http://dx.doi.org/ 10.1103/PhysRevB.111.104509http://dx.doi.org/10.1103/PhysRevB.111.165141http://dx.doi.org/10.1073/pnas.1721495115http://dx.doi.org/10.1103/PhysRevB.99.235149http://dx.doi.org/10.1103/PhysRevB.99.235149http://dx.doi.org/10.1103/PhysRevB.109.085119http://dx.doi.org/ 10.1038/s41586-023-06408-7http://dx.doi.org/ 10.1088/0256-307X/40/11/117302http://dx.doi.org/ 10.1088/0256-307X/40/11/117302http://dx.doi.org/ 10.1103/PhysRevLett.125.257002http://dx.doi.org/ 10.1103/PhysRevLett.125.257002(2020).[49] Hiroyuki Yamase, “Theoretical insights into electronic nematic order, bond-charge orders, andplasmons in cuprate superconductors,” J. Phys. Soc. Jpn. 90, 111011 (2021).[50] M. Hepting, T. D. Boyko, V. Zimmermann, M. Bejas, Y. E. Suyolcu, P. Puphal, R. J. Green,L. Zinni, J. Kim, D. Casa, M. H. Upton, D. Wong, C. Schulz, M. Bartkowiak, K. Habicht,E. Pomjakushina, G. Cristiani, G. Logvenov, M. Minola, H. Yamase, A. Greco, and B. Keimer,“Evolution of plasmon excitations across the phase diagram of the cuprate superconductorLa2−xSrxCuO4,” Phys. Rev. B 107, 214516 (2023).[51] Abhishek Nag, Luciano Zinni, Jaewon Choi, J. Li, Sijia Tu, A. C. Walters, S. Agrestini, S. M.Hayden, Mat́ıas Bejas, Zefeng Lin, H. Yamase, Kui Jin, M. Garćıa-Fernández, J. Fink, AndrésGreco, and Ke-Jin Zhou, “Impact of electron correlations on two-particle charge response inelectron- and hole-doped cuprates,” Phys. Rev. Res. 6, 043184 (2024).[52] Mat́ıas Bejas, Xianxin Wu, Debmalya Chakraborty, Andreas P. Schnyder, and Andrés Greco,“Out-of-plane bond-order phase, superconductivity, and their competition in the t-J∥-J⊥model: Possible implications for bilayer nickelates,” Phys. Rev. B 111, 144514 (2025).[53] P. Prelovšek and P. Horsch, “Electron-energy loss spectra and plasmon resonance in cuprates,”Phys. Rev. B 60, R3735–R3738 (1999).[54] (), in reality, one would expect J⊥/J ∼ 0.1 [65, 66], yielding J⊥/t ∼ 0.03, which, however,may be small enough to be neglected.[55] Mat́ıas Bejas, Andrés Greco, and Hiroyuki Yamase, “Possible charge instabilities in two-dimensional doped Mott insulators,” Phys. Rev. B 86, 224509 (2012).[56] Mat́ıas Bejas, Andrés Greco, and Hiroyuki Yamase, “Strong particle-hole asymmetry of chargeinstabilities in doped Mott insulators,” New J. Phys. 16, 123002 (2014).[57] J. Hubbard, “Electron correlations in narrow energy bands,” Proc. R. Soc. Lond. A 276,238–257 (1963).[58] L. Faddeev and R. Jackiw, “Hamiltonian reduction of unconstrained and constrained systems,”Phys. Rev. Lett. 60, 1692–1694 (1988).[59] K. Sundermeyer, Constrained Dynamics (Berlin: Springer, 1982).[60] Jan Govaerts, “Hamiltonian reduction of first-order actions,” Int. J. Mod. Phys. A 05, 3625–3640 (1990).[61] Adriana Foussats and Andrés Greco, “Large-N expansion based on the Hubbard operator35http://dx.doi.org/ 10.1103/PhysRevLett.125.257002http://dx.doi.org/ 10.1103/PhysRevLett.125.257002http://dx.doi.org/10.1103/PhysRevB.107.214516http://dx.doi.org/ 10.1103/PhysRevResearch.6.043184http://dx.doi.org/ 10.1103/PhysRevB.111.144514http://dx.doi.org/10.1103/PhysRevB.60.R3735http://dx.doi.org/ 10.1098/rspa.1963.0204http://dx.doi.org/ 10.1098/rspa.1963.0204http://dx.doi.org/ 10.1142/S0217751X90001574http://dx.doi.org/ 10.1142/S0217751X90001574path integral representation and its application to the t−J model. II. The case for finite J ,”Phys. Rev. B 70, 205123 (2004).[62] Hiroyuki Yamase, Mat́ıas Bejas, and Andrés Greco, “Electron self-energy from quantumcharge fluctuations in the layered t − J model with long-range Coulomb interaction,” Phys.Rev. B 104, 045141 (2021).[63] Mat́ıas Bejas, Hiroyuki Yamase, and Andrés Greco, “Dual structure in the charge excitationspectrum of electron-doped cuprates,” Phys. Rev. B 96, 214513 (2017).[64] J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, “Evidence forstripe correlations of spins and holes in copper oxide superconductors,” Nature 375, 561–563(1995).[65] S. M. Hayden, G. Aeppli, T. G. Perring, H. A. Mook, and F. Doğan, “High-frequency spinwaves in YBa2Cu3O6.15,” Phys. Rev. B 54, R6905–R6908 (1996).[66] D. Reznik, P. Bourges, H. F. Fong, L. P. Regnault, J. Bossy, C. Vettier, D. L. Milius, I. A.Aksay, and B. Keimer, “Direct observation of optical magnons in YBa2Cu3O6.2,” Phys. Rev.B 53, R14741–R14744 (1996).36http://dx.doi.org/10.1103/PhysRevB.70.205123http://dx.doi.org/10.1103/PhysRevB.104.045141http://dx.doi.org/10.1103/PhysRevB.104.045141http://dx.doi.org/10.1103/PhysRevB.96.214513http://dx.doi.org/ 10.1038/375561a0http://dx.doi.org/ 10.1038/375561a0http://dx.doi.org/10.1103/PhysRevB.54.R6905http://dx.doi.org/10.1103/PhysRevB.53.R14741http://dx.doi.org/10.1103/PhysRevB.53.R14741 Abstract Strong-coupling theory of bilayer plasmon excitations introduction Formalism Results Charge excitation spectra with tz'=0 Overall spectrum Plasmon excitations Charge excitation spectra with tz' =0 Vc dependence of dispersive modes Discussions Conclusions Acknowledgments Large-N formalism of the bilayer t-J-V model Analysis of experimental data with Vc=0 References