Generalized plasma waves in layered superconductors

In a layered and strongly anisotropic superconductor the hybrid modes provided by the propagation of electromagnetic waves in the matter identify two well separate energy scales connected to the large in-plane plasma frequency and to the soft out-of-plane Josephson plasmon. Despite the wide interest in their detection and manipulation by means of different experimental protocols, a unified description of plasma waves valid at arbitrary energy and momentum is still lacking. Here we provide a complete description of generalized plasma waves in a layered superconductors by taking advantage of their connection to the gauge-invariant superconducting phase. We show that the anisotropy of the superfluid response leads to two intertwined hybrid light-matter modes with mixed longitudinal and transverse character, while a purely longitudinal plasmon is only recovered for wavevectors larger than the crossover scale set in by the plasma-frequencies anisotropy. Interestingly, below such scale both modes appears with equal weight in the physical density response. Our results open a promising perspective for plasmonic applications made possible by the next-generation spectroscopic techniques able to combine sub-micron momentum resolution with THz energy resolution.

In a layered and strongly anisotropic superconductor the hybrid modes provided by the propagation of electromagnetic waves in the matter identify two well separate energy scales connected to the large in-plane plasma frequency and to the soft out-of-plane Josephson plasmon. Despite the wide interest in their detection and manipulation by means of different experimental protocols, a unified description of plasma waves valid at arbitrary energy and momentum is still lacking. Here we provide a complete description of generalized plasma waves in a layered superconductors by taking advantage of their connection to the gauge-invariant superconducting phase. We show that the anisotropy of the superfluid response leads to two intertwined hybrid light-matter modes with mixed longitudinal and transverse character, while a purely longitudinal plasmon is only recovered for wavevectors larger than the crossover scale set in by the plasma-frequencies anisotropy. Interestingly, below such scale both modes appears with equal weight in the physical density response. Our results open a promising perspective for plasmonic applications made possible by the nextgeneration spectroscopic techniques able to combine sub-micron momentum resolution with THz energy resolution.

I. INTRODUCTION
Plasmons represent the fundamental excitations of the conduction electrons in metals, and their existence can be easily understood within a classical framework based on Maxwell's equations. Indeed, in a source-free metal bulk plasmons are characterized by zero magnetic field and longitudinal electric field (∇ × E = 0), so that they trivially satisfy Maxwell's equations under the condition of vanishing permittivity 1 . Since the longitudinal electric field couples to density fluctuations the plasma excitation appears also in the charge-density response as collective modes of the electron gas 2,3 . Such a hybrid light-matter mode is also named plasma polariton 4 , especially in the context where spatial confinement at the interfaces between a thin metallic film and a dielectric gives rise to a propagating mode, that can be observed with scanning near-field optical microscopy 5 . When the metal undergoes a superconducting (SC) transitions plasmons characterize also the fluctuation spectrum of the superconducting phase of the complex order parameter formed below the SC temperature T c 6 . Indeed, as originally pointed out by Anderson 7 , the sound-like propagating phase mode of the neutral superfluid 2 is converted into gapped plasma oscillations in a charged superconductor. The appearance of the plasma mode in the spectrum of phase fluctuations is a natural consequence of the fact that the quantum phase of electrons is the variable conjugate to the density [6][7][8][9] . From the theoretically point of view, such an equivalence has been often exploited in the literature to determine the plasma dispersion of a superconductor via the study of the spectrum of the phase mode. The latter can be carried out in a rather elegant and compact way by deriving directly the quantum action for the phase degrees of freedom 6,[8][9][10][11][12][13] , as we will discuss in details in this manuscript.
From the experimental point view there is a strong interest in materials hosting stacks of SC sheets with weak inter-layer SC Josephson coupling 17,18 . For the sake of concreteness, we will model them as xy SC planes stacked along the z direction, see Fig. 1a. The benchmark example of such system is provided by high-temperature cuprate superconductors 26 . The direct consequence of the weak interlayer coupling is that the plasma energy ω z associated with phase (or density) fluctuations among the planes is substantially smaller than the one ω xy connected with in-plane fluctuations. In the SC state an indirect evidence of this is the large anisotropy of the inplane vs out-of-plane penetration depths [27][28][29] , given by λ xy/z = c/ω xy/z . Since the SC gap opening below T c suppresses particle-hole excitations, the low-energy outof-plane Josephson plasmon becomes undamped, and it has been observed long ago via reflectivity measurements as a clear plasma edge emerging below T c at few THz in several cuprates 20,[30][31][32][33] . At finite momentum the dispersion of the Josephson plasmon has been derived by linearizing the sine-Gordon equations for the phase difference among the planes [14][15][16][17][18][19]34 , leading to: where c is the light velocity, γ = ω xy /ω z denotes the plasma-frequency anisotropy, that can be up to two orders of magnitude in some families of cuprates, and k xy , k z denote the in-plane and out-of-plane momentum, respectively. Eq. (1) is sketched in Fig. 1b as a function of k xy at various values of k z . As it is evident from Eq. (1), it represents a regular function of k which tends always to ω z as k → 0. The undamped nature of this propagating mode in the SC state motivated recently several proposals of possible applications. From one side, the inherent non-linearity of the Josephson coupling among CuO 2 planes may enable non-linear photonic applications based on the the use of intense THz laser pulses, as it has been widely addressed both experimentally 18 (1), that is expected to describe the plasma dispersion near the low-energy out-of-plane gap ωz and for momenta much smaller than the scale set by the light cone for the largest plasma frequency. It has been derived [14][15][16][17][18][19] within several papers focusing on the dynamics of the out-of-plane SC phase, in relation to the THz optical response of the Josephson plasmon. Panel (c) shows instead the dispersion (2), that is expected to describe the plasma dispersion near the high-energy in-plane gap ωxy. It is has been derived by studying the generic in-plane and out-ofplane phase or density dynamics in the SC phase 10,11,13,20,21 , and it corresponds to the standard result for the layered electron gas [22][23][24][25] , used to interpret RIXS or EELS experiments, which probes relatively high momenta and energy.
theoretically 17,19,38,39 . From the other side, low-energy Josephson plasmons can be attractive for photonic applications based on near-field nano-optics 4,5,40,41 , that also focus on the regime of few THz and small k ω xy /c momenta.
While optics is interested in the limit of long wavelength, other spectroscopic probes like resonant X-ray scattering (RIXS) and electron energy-loss spectroscopy (EELS) are mostly able to detect the plasma dispersion around ω xy and relatively larger momenta k ω xy /c as compared to the scale set by the light cone 42 . In this regime the effects of anisotropy have been mainly discussed focusing on the emergence of acoustic-like branches of plasma excitations dispersing below ω xy , see Fig. 1c. In this case, the theoretical calculations focusing on the SC state derived the plasma dispersion either looking at the dynamics of the (in-plane and out-of-plane) SC phase [10][11][12][13] or at the dynamics of the density 20,21 . In both cases the standard result is fully analogous to the one found for the metallic state by means of the layered electronic model in the presence of Coulomb interactionsl [22][23][24][25] , and it reduces to a plasma dispersion given by: where an additional ∼ v 2 P k 2 term, relevant to describe the plasmon dispersion with velocity v P at k z = 0, can be added 13 . As one can see in Fig. 1c, at finite k z there is an intermediate regime of k xy values where the plasmon softens below ω xy crossing from a quasi-linear to a ω ∝ √ k behavior, typical of two-dimensional plasmons. So far, such acoustic-like branches have been identified by RIXS in some high-temperature cuprates [43][44][45] , while EELS observed only the high-energy optical branch, with a strong and still not well understood overdamping [46][47][48][49] . Despite the fairly good agreement between Eq. (2) and RIXS experiments, which probe energies around ω xy and large momenta as compared to ω xy /c, there is apparently no connection between Eq. (2) and the expression (1), which should represents the limit of the plasma dispersion when one focuses around ω z and small momenta. Indeed, in contrast to Eq. (1), Eq. (2) is non-analytic as k → 0, since it predicts a continuum of possible values which depend on the angle that the wavevector k forms with z axis. In addition, at k z = 0 and finite k xy , as shown in Fig. 1b,c, Eq. (1) predicts a mode dispersing away from the soft Josephson plasmon ω z while Eq. (2) (with the additional v 2 P k 2 term) predicts a mode around ω xy .
In this paper we provide a general derivation of the plasma dispersion in a layered superconductor which is valid at generic momentum, and we identify a crossover value k c = (ω 2 xy − ω 2 z )/c such that the two previously reported expressions (1) and (2) emerge as the limit of small k k c or large k k c momentum of the full dispersion. To this aim we analyze the problem focusing on the degrees of freedom linked to the SC phase of the order parameter within an effective-action formalism, which naturally implements the spontaneous symmetry breaking below T c . We introduce explicitly the electromagnetic (e.m.) scalar and vector potentials to construct a gauge-invariant quantum action for the coupled system given by the matter (represented by the SC phase) and the e.m. fields. The main technical difference with respect to previous work using the same formalism 10,11,13 , that was able to obtain only the result Eq. (2), relies on the fact that for an anisotropic superconductor the role of e.m. interactions is not exhausted by the Coulomb interaction, which originates from the scalar potential. The reason is that, as one already observes at the level of Maxwell equations [14][15][16]34,38 , in a layered system the e.m. longitudinal and transverse response get intrinsically mixed. The physical reason is that once the superfluid response is anisotropic the induced current is no more parallel to the electric field. Within the context of the SC phase dynamics the anisotropy of the superfluid response leads, in contrast to the isotropic case, to a finite coupling of phase fluctuations to both the longitudinal and transverse components of the e.m. fields. When both couplings are properly included the phase spectrum at arbitrary wave vector carries on the information on a mixed longitudinal-transverse excitation, except for the two special limits (k xy → 0, k z = 0) or (k xy = 0, k z → 0). By introducing an appropriate gauge-invariant phase variables we can derive an analytical expression for the generalized plasma modes valid at generic momentum. In particular, the mode lower in energy, that becomes purely longitudinal only well above k c , if shown to interpolate among the two limiting behaviors provided by Eq. (1) and (2), clarifying then the limits of their validity. In addition, such a formulation allows one to easily extend the result to the case of a non-linear Josephson model, as previously discussed within the context of the out-of-plane Josephson plasmon [14][15][16]34,38 . Finally, in order to make a closer connection with experimental probes sensitive to density correlations 42 , we explicitly derive the densitydensity response in the general case, showing that both generalized plasma modes appear in the spectrum of the physical observables, with an equal weight at the length scale ∼ 1κ c , which is typically of the order of a fraction of micron. At larger momenta k k c the density response closely matches the standard RPA result, and a single quasi-longitudinal plasma mode is visible in the density spectral function. Our results represent a complete gauge-invariant description of generalized plasma excitations in a layered superconductor valid at different energy and momentum scales, providing the benchmark behavior to analyze plasmonic effects in a wide set of experiments ranging from near-field optics to RIXS and EELS.
The plan of the paper is the following. In Sec. II we review the standard derivation of the Gaussian phaseonly action in the isotropic case. The Section starts with an introductory subsection where we outline the theoretical approach used in the manuscript, and a second one where we show its application in the standard isotropic case. Here we show that the plasma mode appears in the different sectors of the matter-e.m. field action depending on the gauge choice, but leading always to the same result. The results of this Section are not new, but their derivation is carried out in a different way with respect to previous work, and it allows us to introduce the formalism relevant for the rest of the paper. In Sec. III we derive the action for the layered 3D case, expressed in terms of the appropriate gaugeinvariant variables. We then show that the longitudinal and transverse e.m. excitations are always coupled, and both excitations appear in the phase spectrum. A simplified description of the problem within the context of Maxwell's equation is also provided in Appendix B. In subsection III C we rephrase the same results in terms of approximate equations of motions, that has been widely used in the recent literature [14][15][16][17][18][19]34 in order to include non-linear effects of the plasmon dynamics. In Sec. IV, we show how the mixed longitudinal-transverse modes found in Sec. III also appears as poles of the properly computed density-density response function. Sec. V contains a general discussion about the results. Further technical details are provided in the Appendixes. Appendix A reviews the main steps leading to the effective action for a superconductor within the functional-integral formalism. Appendix B shows the derivation of the mixed longitudinal-transverse modes from the solution of the Maxwell's equations.

II. EFFECTIVE ACTION FOR THE ISOTROPIC CASE
A. Description of the plasmon via the SC phase Before giving technical details about the derivation of the plasma dispersion for the layered superconductor, it can be useful to briefly outline the general context, starting from the isotropic case. As mentioned in the introduction, the SC transition carries on specific signatures both in the single-particle fermionic excitations spectrum, via the opening of the SC gap below T c , and in the spectrum of the bosonic-like collective excitations. More specifically, while the gap is connected to the equilibrium value of the SC complex order parameter, two modes emerge connected to fluctuations of its amplitude and phase. Since the SC order parameter breaks the continuous gauge symmetry a Goldstone mode is expected, that is directly linked to phase fluctuations of the SC order parameter 6 . There are several ways to see this. At the level of the classical Ginzburg-Landau description this appears as an energetic cost needed to "twist" the phase θ with respect to its equilibrium value. In full analogy with the elastic deformation of a solid body the cost E of a finite phase gradient in space scales with a superfluid "stiffness" D s , such that In the simplest Galilean-invariant case D s can be expressed as the ratio between the density of superfluid electrons and their mass, D s = n s /m, with n s = n at T = 0. As the length scale of spatial phase deformations goes to zero E vanishes. This can be interpreted as a mode at zero energy for k = 0, that is what expected for a Goldstone mode. On the other hand, to really define a propagating mode, i.e. to establish a frequencymomentum dispersion, one needs to include dynamical effects, not present in the classical Ginzburg-Landau description. A very elegant and powerful technique relies 6 on the explicit construction of a quantum analogous of the Ginzburg-Landau model by starting from a microscopic interacting model for electrons. The basic idea is to start from a fermionic model with a BCS-like interaction term and decouple it via the Hubbard-Stratonovich procedure by introducing two effective bosonic fields which play the role of the order-parameter amplitude and phase. By further integrating out explicitly the fermions one obtains a quantum model expressed in terms of the collective SC variables, that can be expanded in principle up to arbitrary powers in the collective-mode fields, linking the phenomenological couplings of the Ginzburg-Landau model to fermionic susceptibilities, which contain all the information on the microscopic physics under investigation. Such a procedure is outlined in Appendix A, and further details can be found in several books and papers on the subject, see e.g. Ref.s 6,[8][9][10][11][12][13] . By retaining only Gaussian terms in the fluctuations one defines the spectrum of the collective modes, equivalent to compute the response functions at RPA level in the usual diagrammatic approach to fermionic models 2,6 . Within this quantum generalization of the classical Ginzburg-Landau theory one finds that the density appears as conjugate variable of the phase 6,8,9 . As a consequence, within the quantum phase-only action the energetic cost to perform a phase gradient in time is controlled by the charge compressibility κ 0 8,9 , and Eq. (3) is replaced at T = 0 by the quantum action: For weakly-interacting neutral systems κ 0 in the static long-wavelength limit can be approximated with the density of states at the Fermi level, and in Eq. (4) one recognizes the so-called 7 Anderson-Bogoliubov sound mode: where v 2 s = D s /κ 0 is the sound velocity. The identification of the density as conjugate variable of the phase is a general result, which holds both for neutral superfluids and for superconductors, and has profound implications for the latter. Indeed, in the standard description of charged systems, density-density interactions are expected to be mediated by the longrange Coulomb potential. In the literature based on the effective-action formalism the difference between neutral and charged superfluids is then usually encoded 6,9-13 via an additional Hubbard-Stratonovich field which decouples density-density fermionic interactions, so that its integration dresses the charge compressibility at RPA level by the Coulomb potential 7 . The term κ 0 in Eq. (4) is then replaced in Fourier space by κ(k) = κ 0 /(1 + V (k)κ 0 ), where V (k) is the Coulomb potential in generic D dimensions. Since for k → 0 one has κ → 1/V (k) the spectrum of the phase mode, that reflects the one of density fluctuations, identifies now a plasma mode, whose energy vs momentum dispersion depends on the dimensions. In the standard isotropic three-dimensional (3D) case one recovers the well-known dispersion where ω 2 p = 4πe 2 D s is the isotropic plasma frequency. The procedure outlined above explains why in the SC state the spectrum of the collective plasma excitations can be derived in a relatively simple way by analyzing the SC phase fluctuations. However, this analysis is not enough if one wants to determine the spectral weight of the plasma mode, that can only be determined by direct derivation of the density-density response 6,10,12,13 , which couples to external probes. Indeed, on a more general ground, the SC phase itself does not represent a physically observable quantity. Once the matter is coupled to the e.m. field, the physical observables are related to the gauge-invariant four-vector ψ µ , whose temporal and spatial components read: where φ and A are the scalar and vector potentials respectively. In contrast to the phase alone, the observables (7) are invariant under the simultaneous gauge transformation of the e.m. potentials (that leaves the physical electrical and magnetic field unchanged) and of the SC phase: In this picture one realizes that analyzing the spectrum of the phase mode in the superconductor is completely equivalent to solve the problem of the e.m. wave propagation. More specifically, the frequency-momentum relation (6) corresponds to the one of the longitudinal component of the electric field, which couples by Maxwell's equations to charge fluctuations and then to phase fluctuations in a superconductor. In general, e.m. waves can also be transverse, but in the isotropic system they are not coupled to phase fluctuations. Once again, this result can be easily understood already at classical level. In isotropic and homogeneous systems the minimal-coupling substitution (7) into the gradient term (3) leads to The explicit coupling among the phase and the gauge field in Eq. (4) is the term dxD s (∇θ) · A, that for a constant D s can be written after integrating by part as dx D s θ(∇ · A). As a consequence the phase couples only to the longitudinal component of the e.m. field, and one does not see any signature of the transverse e.m. mode in the phase spectrum. Such a decoupling of the longitudinal and transverse sectors makes the description of the plasma mode completely equivalent in the two approaches: via the matter, where the charged nature of the system enters via the Coulomb interaction for the density, or via e.m. potentials, where one derives directly the equations of motions for the e.m. fields.
The situation becomes radically different for anisotropic systems, and in particular for layered superconductors. In this case, which is relevant for several classes of unconventional superconductors, as e.g. cuprates or pnictides, the phase stiffness in the plane D and in the direction perpendicular to the planes D ⊥ is anisotropic, with D ⊥ D [27][28][29] . Within the language of Eq. (9) above, this implies that the phase can also be coupled to the transverse vector potential. This is the central physical effect which motivates the present work, aimed to explain what are the consequences of such a mixing of the transverse and longitudinal modes for the spectrum of plasma waves. Indeed, as we mentioned in the introduction, so far the effect of this mixing has been included in an approximated way in the work [14][15][16] aimed to derive the dispersion (1) of the soft Josephson plasmon, but it is completely missing in any approach extending the RPA to the layered electron gas [22][23][24][25] or to the layered superconductors 10,11,13,20,21 . Indeed, in the latter case, which leads to the expression (2) for the layered plasmon, one takes into account only the effect of the Coulomb potential, which is linked to the longitudinal component of the e.m. field.
In Sec. (III) we will explicitly show how the transverselongitudinal mixing appears in the quantum phase-only model for the layered superconductor, and how one can derive from it an analytical expression for the generalized plasma dispersion which interpolates between the two limits (1) and (2), clarifying at the same time their range of validity. To get a better insight into the link among the phase mode and the plasma wave we will first review in the next subsection the isotropic case. In contrast to previous work 6, [9][10][11]13 where the e.m. field only enters via the Coulomb density-density interaction, we introduce here the coupling (7) of the SC phase to the internal e.m. field, that we further integrate out. We then show how the same plasmon dispersion is obtained regardless the gauge choice for the e.m. field, as expected. This approach will then be extended to the layered case in Sec. (III).

B. Formulation of the problem with the internal e.m/ fields
As a starting point we use the imaginary-time equivalent of the quantum action (4) for an isotropic 3D superconductor, rewritten in Fourier space. As mentioned above, its derivation relies on a rather standard technique 6,[8][9][10][11][12][13] , whose details are given in Appendix A. We then have: where q ≡ (iΩ m , k) is the imaginary-time 4-momentum, Ω m = 2πmT are bosonic Matsubara frequencies, κ 0 is the bare compressibility, and D s is the superfluid stiffness. In the following we put = k B = 1. In ordinary Maxwell equations the source for the e.m. fields are provided both by external charge/currents and by the induced internal charge/current fluctuations in the matter. In the superconductor the latter will be described by the SC phase, coupled to the internal e.m. fields via the minimal coupling (7). To account for the e.m. fields we then need to include first the free contribution, which reads 6 : Here we introduced the longitudinal , where the electric field is expressed in the imaginary-time formalism as It is worth noting that in order to have a definition of E analogous to the one valid for real time one should assume that φ is purely imaginary, i.e. one should replace φ → iφ. In this case, by defining the imaginary-time electric field as E ≡ − 1 c ∂ τ A − ∇φ the action for the free e.m. field would read ε|E| 2 +|B| 2 8π , which is the usual expression for the energy density. Such a rescaling of the scalar potential would also make the quadratic term in the scalar potential arising from (∇φ) 2 positive defined, as required to perform the Gaussian integration. To make notation more compact we will not explicitly rescale the potential in what follows, but we will implicitly assume that a formal definition of the Gaussian integration in the imaginary-time formalism requires such a regularization. In order to include the ionic screening, we also introduced the background dielectric constant ε. As mentioned above, in the superconductor 6,9,10,12 the coupling of the e.m radiation with the matter can be easily implemented by using the SC phase, via the minimal coupling substitution (7), which reads in momentum space and Matsubara frequency: By replacing Eqs. (12) into the action (10) one obtains a light-matter coupling term S θAµ [θ, A µ ] and a renormalization of the bare e.m. action, so the total action reads: where S e.m. [A µ ] describes the long-wavelength propagation of light through matter. For the isotropic system this term reads: In the effective-action formalism the frequencymomentum dispersion obtained as solution of the Maxwell's equations appears as the zero of the determinant of the matrix associated with the coefficients of the fields in the Guassian action, once the analytical continuation iΩ n → ω + iδ is performed. As a consequence, one clearly sees that the inclusion of the matter has two well-known effects on the e.m. fields: (i) the scalar potential displays the usual Thomas-Fermi screening, where k 2 T F = 4πe 2 κ 0 /ε; (ii) the propagating transverse mode is gapped out 50 at the isotropic (screened) SC plasma frequency and the light velocity is renormalized asc = c/ √ ε. In the typical field-theory language, the Anderson-Higgs mechanism 50 manifest as a "mass" term ω 2 p A 2 T in the e.m. action (14), that, in the static limit, only survives below T c , where D s = 0. Finally, the coupling between the phase and the e.m. potential reads: The last term of Eq. (17) shows explicitly that in the isotropic case phase fluctuations only couple to the longitudinal component of the vector potential A L , as we mentioned above. In addition, the use of Eq. (12) clearly guarantees that the total action is invariant under the gauge transformation (8), which reads explicitly in Matsubara formalism: λ(q) being an arbitrary function of the Fourier momenta. Such a gauge freedom also implies that only the gauge-invariant combination (12) represents a physicallyobservable quantity describing matter properties, while the information carried out by the SC phase alone acquires a different meaning within different gauge choices. This means that even though the phase-only propagator is always linked to the plasma mode, it will appear in different forms depending on the gauge. To clarify this point, let us start from Eq. (13) and let us derive the action for the phase degrees of freedom by integrating out the e.m. potentials. A first natural gauge choice is the so-called Coulomb gauge ∇ · A = 0, i.e. A L = 0. Indeed, as noticed below Eq. (17), in the isotropic case the phase fluctuations only couple to A L , so that in the Coulomb gauge only the coupling between the phase and the scalar potential survives. By integrating out φ one is then left with: where we defined α as where λ D is the Debye screening length. The result (19) is formally identical to the one widely discussed in the previous literature 6,9,10,12,13 , and obtained by adding a Coulomb-mediated density-density interaction term in the starting microscopic fermionic model. This is a natural consequence of the fact that in the Coulomb gauge only the coupling of the phase to the scalar potential is relevant.
An alternative but completely equivalent approach can be followed integrating out the e.m. fields in the Weyl gauge, where φ = 0. In this case only the coupling to A L survives in Eq. (17) and one obtains: We immediately see that both (19) and (21) identify a collective excitation as a pole of the θ(q) fluctuations, which are controlled by the inverse of the |θ(q)| 2 term. After analytical continuations to real frequencies this is given by: that is the usual dispersion of the longitudinal plasma mode in a superconductor. As usual, the plasmon velocity c p = √ αω p =c(λ D /λ), where λ is the London penetration depth, is much smaller than the light velocityc of the transverse wave in the medium, since λ D /λ 1 even for unconventional superconductors. It is worth noting that the velocity of the plasmon in Eq. (22) is not the same as the one obtained for the normal metal 3 , that would correspond to ω 2 = ω 2 p (1 + (9/5)α|k| 2 ). This is due to the fact that to correctly account for the plasma dispersion one should also account for the |k| 2 corrections to the BCS density-density response function in the SC state, that appears as a coefficient of the (∂ t θ) 2 term in the phase-only action, while in Eq. (4) only its long wave-length limit κ 0 has been retained. To simplify the notation we shall not perform explicitly such an expansion, that is anyway irrelevant to the physical effects we want to discuss in the present manuscript. Indeed, it only gives a small quantitative difference that can be included if one is interested into a detailed quantitative comparison with the experiments. While in Eq. (19) the plasma mode appears in the spectrum of ∇θ fluctuations, in the Weyl gauge it appears in the spectrum of ∂ τ θ fluctuations. This is a direct consequence of the fact that the phase variable does not represent a true physical observable. In addition, it describes, as expected, only the longitudinal mode. For both gauge choices the transverse e.m. mode is described by the two remaining degrees of freedom of the e.m. field in Eq. (13). A very elegant and convenient way to derive simultaneously the energy-momentum dispersion for all the transverse and longitudinal e.m. modes relies on the use of the gauge-invariant physical observables (7). For example one can use the spatial component of the gauge-invariant (g.i.) phase difference: and set to zero the scalar component via a proper gauge fixing. The equivalent of Eq. (13) in the new variables reads: By introducing again the longitudinal ψ L = (k · ψ)k and the transverse ψ T = ψ−ψ L = (k×ψ)×k components we see that only ψ L couples to the phase θ. After integrating it out one then finds: From Eq. (25) one immediately sees that the three components of ψ describes all e.m. modes, as given by the poles of the longitudinal and transverse propagators. Indeed, after analytical continuation, one finds the solution (22) in the longitudinal sector and the solution (15) in the transverse sector, respectively. Such a result suggests that in the anisotropic case where longitudinal and transverse modes get mixed a description in terms of the physical fields (23) can be more convenient, as we shall see indeed in the next section. In addition, this choice will make the extension of the Gaussian action to a nonlinear Josephson model straightforward, since it will only affect the mass term for the ψ field, as we will see in details in Sec. III C.

III. EFFECTIVE ACTION FOR A LAYERED 3D SYSTEM
ol(x,y)  Fig. 1c.
A. Derivation of the plasma dispersion for the anisotropic system To discuss the layered case we should consider an extension of Eq. (10) where the stiffness is anisotropic. By using the notation of Fig. 1a we will denote with D xy and D z the in-plane and out-of-plane stiffness, respectively. Eq. (10) can be straightforwardly generalized 10,11,13 to: where k xy ≡ k 2 x + k 2 y is the in-plane momentum. Notice that, despite of the discrete out-of-plane structure of the system, we still use a continuum anisotropic model, as valid when k z 1/d. Indeed, the breakdown of the standard RPA approximation and the mixing between longitudinal and transverse degrees of freedom, which will be the main topics this section, occur at very small momenta k, which are suitably accounted for in the continuum limit. The anisotropy of the stiffness has important consequences when the electrons are coupled to the e.m. fields via the minimal-coupling substitution (12). Indeed, the action (14) describing the e.m. field in the matter gets modified as: where we defined the two plasma frequencies: In principle, in Eq. (27) we could also consider an anisotropic background dielectric constant, in order to take into account the different ionic screening occurring along the in-plane and the out-of-plane directions. The e.m. energy density would then read −εα|Eα| 2 +|B| 2 8π , with ε α = ε xy for α = x, y, and ε α = ε z for α = z. However, this effect does not modify substantially the physics of the system, therefore for the sake of simplicity we take ε isotropic. Finally, the coupling of the e.m. fields with the phase variable is an anisotropic version of Eq. (17), that is now expressed as: As one can see, the anisotropy of the SC phase stiffness has two main consequences in the description of the e.m. response: (i) the massive terms in the vector potential of Eq. (27) become anisotropic; (ii) the inclusion of the matter via the SC phase shows that there is no way to decouple completely the longitudinal and transverse modes at arbitrary wave vector k. The latter result emerges more clearly by close inspection of Eq. (29). While in the isotropic case the equivalent term D s k · A of Eq. (17) implies that θ only couples to the longitudinal component of the vector potential, in the anisotropic case the combination D xy k xy · A xy + D z k z · A z is different from zero even in the Coulomb gauge where k · A = 0. This is the crucial point that has been missed in previous work where the effect of the anisotropy has been only included in the phase stiffness, but not on the coupling to the gauge potential. At a more general level, the physical effect we are implementing here is the fact that the superfluid current is no more parallel to to electric field, so even a purely longitudinal field always induces a transverse current and viceversa. Such an interpretation is straightforward when the problem is studied via Maxwell's equations only, as we discuss in Appendix B.
To better understand the consequences of the coupling between transverse and longitudinal modes, let us first recall the derivation of the standard result (2). If one uses the Coulomb gauge and retains in Eq. (29) only the coupling to the scalar potential φ, the effect of its integration is a straightforward generalization of Eq. (21), provided that one accounts for the anisotropy of the stiffness of Eq. (26). The action for the phase would then read: The plasmon dispersion obtained after analytical continuation from the phase propagator in Eq. (30) in the limit α → 0 is exactly the result (2), i.e.
where η denotes the angle between the k vector and the z axis. Eq. (31) is plotted in Fig. 2 for the values ω xy = 1 eV and ω z = 0.05 eV, which lead to an anisotropy γ = ω xy /ω z comparable to typical values in cuprates. The result (31) is analogous to the RPA one derived by including only the Coulomb potential. This approach, that we will refer to as standard RPA in what follows, has been adopted in previous work both in the SC [10][11][12][13]20,21 and in the normal 22-25 state. Within this scheme it is also possible to generalize the expression (31) by retaining the lattice periodicity in the z direction 10,13,22-25 . As noticed in the introduction, the expression (31) has a singular limit at k = 0, since its value depends on the angle η. As we shall see below, such a singularity is removed by inclusion of the coupling of the phase to the vector potential A. Finally, for the sake of the following discussion let us also write explicitly the transverse mode obtained in the same approximation, where one neglects the coupling between A T and θ. In this case the dispersion of the transverse e.m. wave is simply encoded in the A 2 T term of Eq. (27) and it reads: To fully account for the coupling of the SC phase to both A L and A T let us take advantage of the results of Sec. II and let us introduce the gauge-invariant phase variables (23). To simplify the notation, we can assume without lack of generality that the in-plane momentum is along the x direction. In this case the ψ y component decouples from the phase in Eq. (29): it describes a pure massive in-plane transverse mode. The remaining two components are coupled, and with lengthy but straightforward calculations one can derive the generalization of Eq. (25): We can now compute the propagators ψ α (q)ψ β (−q) , whose poles identify the collective-mode excitations. From |ψ y (q)| 2 one immediately finds the dispersion of the massive in-plane transverse mode ω 2 = ω 2 xy +c 2 |k| 2 , that simply reflects the standard isotropic-case result (15). The fluctuations of ψ x and ψ z are coupled, so the collective modes are given by the determinant of their 2 × 2 matrix, once the analytic continuation iΩ m → ω +i0 + has been performed. The dispersion of the collective excitations is then obtained as solutions of a quartic characteristic equation, which reads, in the α 0 limit: The solutions of Eq. (34) are: Eq. (35) is the first central result of this paper: it provides the general dispersion of the e.m. excitations of a layered system for any momenta. The two solutions (35) are shown in Fig. 3 as two-dimensional maps, while in Fig. 4 we report the modes as a function of |k| for selected value of the angle η that k forms with the z axis. First of all, we notice that ω + and ω − are regular functions as k → 0. Indeed, assuming that ω xy > ω z , as it  happens typically in layered materials where planes are weakly coupled, one immediately sees that: regardless of the direction along which such limit is taken, as it is better shown in Fig. 4. At finite k both solutions give an increasing value of the plasma excitation, that is limited above by ω xy for the ω − solution, while it increases rapidly for the ω + solution. Such a behavior can be better understood in Fig. 4 where we plot for comparison also the standard RPA results ω 2 L and ω 2 T derived neglecting the coupling of the phase to the transverse field, see Eq. (31) and (32) respectively. Indeed, as |k| increases the two solutions (35) reduce to:  By closer inspection of Eq. (35) one sees that the crossover among the two regimes occurs around a critical value k c of the momentum set by the ratio between the plasma-frequency anisotropy and the light velocity: Indeed, as soon as k k c the square-root term in Eq. (35) can be expanded in powers of the small parameter k c /k and one easily recovers the two analytical expressions of the purely longitudinal and transverse modes ω L and ω T , respectively. To get an idea on the order of magnitude of the momentum scale (38), one should consider that c 0.19 eVµm. Thus, considering that in layered materials as e.g. cuprates it is usually ω z ω xy , with ω xy 1 eV, one sees that as soon as k k c ∼ 5 µm −1 the standard RPA result (31) is recovered, as shown in Fig. 4. As a consequence, for experiments like EELS or RIXS, which measure the plasma dispersion at moment of the order of 1/a ∼ 0.1Å −1 , with a lattice spacing, and energies of the order of the eV, the mixing between longitudinal and transverse modes is not quantitatively relevant. Nonetheless, the discrepancy between ω − and ω L becomes crucial in order to understand the radically different description of the low-energy plasma mode provided within the context of non-linear Josephson plasmonic 17,18 . These experiments are carried out with THz radiation approximately resonant with the lowenergy mode ω ω z and which, at the boundary with the medium, is polarized along the z axis. If one considers the solution of Eq. (34) for ω ω z , one can approximate ω 2 − ω 2 xy −ω 2 xy since ω xy ω z . Physically, this is equivalent to neglect the (∂ψ x /∂τ ) 2 term in the first line of Eq. (33) with respect to the ω 2 xy ψ 2 x term, i.e. to assume a stationary in-plane current, that can be a reasonable approximation at the frequency scale of the soft Josephson plasmon. From the point of view of the eigenomdes of Eq. (33), this approximation turns the quartic equation (34) into a quadratic one, that can be easily solved leading to Eq. (1) mentioned in the introduction, that we can also rewrite as: where we introduced the in-plane and out-of-plane penetration depths: As shown in Fig. 4, Eq. (39) accounts indeed for the correct behaviour of ω − at energies around ω z and small momenta. It should be noticed that for optical THz probes the relevant value of the momentum k in Eq. (39) is of the order of 10 to 100 cm −1 , so well below the threshold (38), i.e. in the region where the standard RPA approximation fails. This explains the apparent disagreement between the two expressions usually quoted in the literature for the same problem, studied within the context of different experimental probes.

B. Longitudinal-transverse mixing at generic wavevector
To gain further insight in the mechanism connecting the general solutions ω ± to the results obtained with the standard RPA approximation it is instructive to have a closer look to their polarization with respect to k, in order to understand how the longitudinal-transverse mixing changes around k c . Let us first analyze the lowmomentum limit k → 0. In this regime the eigenvectors corresponding to the two solutions ω ± read: As one clearly sees, in the k → 0 limit the ψ − solution is always polarized along z, while ψ + is polarized along x, see Fig. 5a, explaining why the plasma modes always reduce to the ω z/xy values, without the continuum of k → 0 values predicted by the standard RPA, see Fig. 2. On the other hand, this also implies that at small momenta the eigenmodes have a mixture of longitudinal and transverse character. This can be better seen by rewriting the two solutions at arbitrary value of the momentum by using explicitly the longitudinal ψ L and transverse ψ T components of the g.i. phase: (42) Since we are setting k in the xz plane, one of the transverse components ψ xz T = (ψ T ×ŷ)×ŷ lies in the xz plane along the directionk ×ŷ, while the other ψ y T = (ψ T ·ŷ)ŷ is perpendicular to it, see Fig. 5. As seen in Eq. (33), ψ y is the polarization of the pure transverse mode not involved in the mixing. The orthogonality between ψ − and ψ + is ensured by the fact that ω 2 . Eq. (42) highlights how, at generic value of the momentum, ω + and ω − describe two modes with no pure (longitudinal or transverse) character. On the other hand as soon as k k c and the modes reach the standard RPA value, see Eq. (37), ψ − becomes a purely longitudinal mode while ψ + is a purely transverse one, see Fig. 5b, consistently with the expectation that the standard RPA approximation leads to pure longitudinal/transverse modes. A remarkable exception to the mixing is provided by the case of propagation completely in plane (k z = 0) or completely out-of-plane (k x = 0), since in this case the ψ x ψ z coupling of Eq. (33) vanishes and one obtains purely longitudinal/transverse modes. For instance when k z = 0 the momentum is along x, so ψ + , which is also aligned along x (see Eq. (41)), describes a purely longitudinal mode approaching ω xy as k x → 0. Conversely, when k x = 0 the momentum is along z, so that in this case the longitudinal mode is represented by ψ − , and its dispersion approaches ω z as k z → 0.

C. Comparison with the sine-Gordon equations of motion
As mentioned in the introduction, the approximated expression (39) for the small k limit of ω − has been widely used in the recent literature [17][18][19] to study the dynamics of the Josephson plasmon in the presence of non-linear effects due to intense THz fields. Within this context, it has been derived the equation of motion for the gauge-invariant variable ψ z of Eq. (7). The basic equation reads: where absence of dissipation is assumed, for the sake of simplicity. The explicit appearance of the sin ψ z term makes the equation non-linear. As mentioned above, this result shows how the use of the g.i. phase variable repre-sents a powerful and elegant way to obtain a straightforward extension of the Gaussian model to the non-linear Josephson model, that represents a crucial step in order to describe non-linear optical effects. On the other hand, the linearized version of Eq. (43), obtained when sin ψ z ψ z , admits a wave-like solution ψ z ∝ e i(ωt−k·r) such that the frequency ω and the momentum k satisfy the dispersion relation (39). As a consequence, in the linear regime Eq. (43) provides again an approximation for the ω − solution at low energy and momenta. It is then interesting to understand how such a non-linear extension to the Josephson model can be obtained starting from the more general formalism developed so far.
A simple way to see the analogy is to rewrite the general result (33) in real space. Consistently with the derivation of Eq. (35), we shall consider here the α = 0 limit. In addition, since Eq. (43) is written in real time we will convert Matsubara frequencies to real frequencies via analytical continuation and we will replace dτ = idt. The resulting real-time actionS ani then reads: Notice that we promoted the ω 2 z ψ 2 z term of Eq. (33) to a cosine-like term −2ω 2 z cos ψ z , so that interacting terms for the phase are included beyond the Gaussian order. As usual, this is physically motivated by the idea that a Josephson-like coupling set by the out-of-plane phase stiffness D z exists for the phase in neighboring layers, leading to an effective XY model for the phase degrees of freedom. The main consequence of the presence of a cosine term is that the resulting phase-only model admits naturally a non-linear current along the z direction, I z ∝ ω z sin ψ z , that is crucial to account for the nonlinear optical response measured at strong THz fields aligned along the z direction in cuprates 35,36 . As discussed in Ref. 12 , the interacting terms for the phase derived microscopically can differ from the one obtained within the simple XY model. Nonetheless, the XY model provides a reasonable starting point to account for non-linear effects in the z direction, as discussed recently in Ref. 39 . Finally, to account for the layered structure the gauge-invariant phase variable ψ should be promoted to a layer-dependent variable ψ(r, z = nd) → ψ n (r, t), with r = (x, y), n layer index and d spacing between layers. The final result is completely analogous to Eq. (44), provided that one discretizes the integration along z as dx → d n dr and defines a discrete version of the derivative along the z direction, so that ∂ z f ≡ 1 d ∂ n f = fn+1−fn d . We will then retain for simplicity the continuous notation in what follows.
Once rewritten the action in real space we can obtain the Euler-Lagrange equations of motion by functional derivatives with respect to ψ x and ψ z . By simple algebra we then obtain two coupled equations which describe the dynamics of the mixed longitudinal-transverse modes: As one can easily check, when sin ψ z ψ z one recovers two coupled linear equations which can be solved with wave-like solutions propagating with frequency ω and momenta k satisfying the same Eq. (34) derived above. Thus in the linear regime, as expected, one recovers the same expression (35) for the ω ± collective modes derived by the Gaussian action (33). On the other hand, if one is interested to study the collective dynamics of ψ z for frequencies around ω z one can get an approximate equation of motion by noticing that for ω ω z ω xy the first time-derivative term of Eq. (45) is of order (ω z /ω xy ) 2 as compared to the second one, and can then be neglected. Thus one simply deduces from Eq. (45) that: where we introduced the penetration depths (40). By applying the (1 − λ 2 xy ∂ 2 z ) operator to Eq. (46) and using Eq. (47) one obtains exactly Eq. (43). In this way, we reconciled all the expression discussed in the previous literature, by clarifying also the limit of validity of the various approximations.
A last comment is in order about the role of the α|k| 2 terms, neglected while deriving the ω ± expressions in Eq. (35). In the isotropic case, these terms are crucial in order to get the plasmon dispersion, see Eq. (22). Since α = λ 2 D coincides with the Debye length squared, see Eq. (20), its effect is negligible as compared to the much larger variations in k between ω z and ω xy already described by Eq. (35). On the other hand, including these corrections into the general solutions ω ± is a matter of simple algebra, and in the standard RPA regime one would just recover an additional α|k| 2 dispersion into Eq. (31).

IV. GAUGE-INVARIANT DENSITY-DENSITY RESPONSE INCLUDING RELATIVISTIC CORRECTIONS
So far, we investigated the nature of the mixed longitudinal-transverse mode by looking at the dynamics of the gauge-invariant phase variable, and we showed a failure of the standard RPA approach to capture the k → 0 limit of the plasma dispersion. It is then worth investigating how such generalized plasma modes appear in the physical gauge-invariant density-density response function, which carries information on the longitudinal degrees of freedom of the system.
To clarify the procedure, we first outline the derivation of the density-density response function for an isotropic SC system. The simplest way is to add an auxiliary scalar field δφ coupled to the density operator into the microscopic hamiltonian, so that the density-density response can be obtained by functional derivative of the total action with respect to δφ. In practice, the idea is to derive an action quadratic in |δφ| 2 after integration out of all the other variables, so that the density-density response function K will be the coefficient of the |δφ| 2 term in the final action, i.e.
where S[δφ] has been obtained by integrating out all other degrees of freedom except then the auxiliary field.
In the isotropic case the result is straightforward. Indeed, once known the phase-only action (10) the scalar perturbation δφ enters the Gaussian phase-only action via the usual minimal coupling substitution (12). One is then left with an effective gaussian model for the variables θ, δφ and φ, where the action S δφ adds to the contribution (24) computed before. The action S δφ contains a bare quadratic term, which accounts for the bare compressibility, and the coupling of the phase with the scalar perturbation. These read: As we did before, we can introduce the gauge-invariant phase variable ψ and integrate out explicitly θ. This leads to a dressing of the bare compressibility term which multiplies |δφ(q)| 2 in Eq. (49) and generates a coupling among δφ and ψ: . (50) As expected, the scalar perturbation couples only to the longitudinal component ψ L of the gauge-invariant phase difference, carrying on the information on the longitudinal modes. Once ψ L is integrated out, we get the fully dressed density-density response function from Eq. (48). Neglecting the plasmon dispersion (i.e. α 0) and performing the analytic continuation iΩ m → ω + iδ it reads: Eq. (51) displays a singularity at the plasma frequency ω p . One can easily check that Eq. (51) is perfectly equivalent to the standard RPA result K RP A (q) = χ0(q) 1+V C (k)χ0(q) obtained in the normal state 3,6 , where χ 0 (q) is the Lindhart function for the metal. Indeed, in the limit ω/k 1 the Lindhart function can be approximated as χ 0 (q) − n m |k| 2 ω 2 , that replaced into the standard RPA formula gives exactly Eq. (51). Since in this frequency/momentum regime one would expect that the SC susceptibility has the same behavior of the normal-state one, one can conclude that Eq. (51) is equivalent to the standard RPA result for the isotropic superconductor. By retaining terms in α|k| 2 in Eq. (50) one could get also the plasmon dispersion. However, the extension of Eq. (51) will not give the correct damping of the plasmon, since it has been derived taking directly the static limit of the density response, i.e. κ 0 = χ 0 (ω = 0, k → 0). To get the full plasmon spectral function one should retain the full frequency and momentum dependence of the bare density-density response χ 0 , in order to recover a plasmon damping when ω p (k) overlaps the particle-hole continuum 3,6 . In order to generalize the previous derivation to the anisotropic case we should introduce, in full analogy with Section III, the in-plane and out-of-plane superfluid stiffness D xy and D z . This leaves unchanged the coupling between δφ and ψ L in Eq. (50), but it modifies the term in |ψ L | 2 . Once again, the standard RPA result does not account for this effect and it only introduces the anisotropy of the bare Lindhart function. For an anisotropic metal this would be equivalent to approximate χ 0 (q) − n ω 2 where the anisotropic in-plane m xy and out-of-plane m z masses identify the anisotropic plasma frequencies ω 2 xy/z = 4πe 2 n/m xy/z . As a consequence, the standard RPA result for the anisotropic case is fully equivalent to Eq. (51), provided that ω 2 p is replaced by its layered version ω 2 L defined in Eq. (31), i.e.
As discussed before, we expect that such an approximation fails in the k k c regime. Indeed, by extending the procedure outlined above for the isotropic case and by integrating out the longitudinal g.i. component ψ L at arbitrary momentum we obtain the following generalization of Eq. (51) in the k k T F limit: From Eq. (53) one sees that the poles of K (ani) are the generalized plasma modes ω + and ω − at generic k. This is consistent with the fact that both have a finite longitudinal projection. Also in this case the full computation of the density spectral function requires a careful evaluation of its imaginary part, that is beyond the scope of the present manuscript. Nonetheless, the result (53) already provides us with the direct access to the total spectral weight I ± (k) allocated in each ω ± (k) pole as a function of the momentum k, as given by the imaginary part of Eq. (53). These can be readily obtained as: The momentum dependence of I ± (k) is shown in Fig. 6 for two values of the angle η between k and the z axis, corresponding to the panels (a) and (c) of Fig. 4. Notice that in order to compare the relative spectral weights at a fixed k we did not include in Eq. (54) the overall |k| 2 factor present in the density-density response (53). As a consequence, in the k k c regime where the ω ± solutions are different from the standard RPA ones ω T /L , the upper pole at ω + has always a larger spectral weight than the lower one ω − , due to the overall ω 2 factor in Eq. (53). Conversely, in the k k c limit, where ω − reduces to the pure RPA longitudinal mode and ω + to the pure RPA transverse mode (see Eq. (37)) only the lower pole at ω − survives. Indeed, in this regime one sees from Eq. (54) that the factor ω 2 + − ω 2 T in the numerator of I + vanishes, while the factor ω 2 − − ω 2 + ω 2 − − ω 2 T at the denominator of I − cancels out exactly the weighting factor ω 2 − − ω 2 T in the numerator, so that I − approaches the spectral weight of the standard RPA longitudinal mode, which is defined from Eq. (52) as I RP A = ω L (k)/2. In other words, in this regime the density-density response approaches K ani RP A in Eq. (52), obtained in the standard RPA approach. Conversely, at k k c Eq. (53) represents a rather unexpected result, i.e. the simultaneous signature in the density response of two generalized plasma modes, located at two well separated energy scales, with the initial predominance of the higher-energy solution ω + .

V. CONCLUSIONS
In the present manuscript we provided the full description of the generalized plasma modes in bulk layered superconductors, filling the knowledge gap among previous results, that focused on specific regions of the energy/momentum spectrum of the plasmons. The main physical mechanism behind our findings is the anisotropy of the superfluid response in layered superconductors, that leads to an induced superfluid current which is no more parallel to the electric field. As it is already evident at the level of classical Maxwell's equations (Appendix B), this induces a mixing between transverse and longitudinal components of the e.m. fields, usually absent in isotropic systems. To rephrase the same effect in an other way, in the anisotropic system density fluctuations are not only generated by a scalar potential but also by a vector potential, and within the context of the superconductor it appears as a finite coupling of the SC phase to the transverse vector potential. Such effects becomes typically relevant for wavevectors that, in strongly anisotropic systems (where i.e. ω xy ω z ), are smaller than k c ∼ ω xy /c. Since this scales vanishes as c → ∞ these effects appear as relativistic corrections. At momenta well above the crossover scale the modes acquire an almost pure longitudinal or transverse character, and the anysotropy only manifests in the emergence of acoustic-like plasmon branches reminiscent of the plasma modes in purely 2D systems.
The above results have been derived by taking advantage of the fact that in a superconductor the plasma modes can be conveniently studied by means of the dynamics of the SC phase of the order parameter. However, previous approaches focused on apparently different models, making it difficult to establish a link between the different expressions for the plasma dispersion discussed in the literature. On this respect, the aim of the present work is twofold. From one side, we provided an analytical expression for the generalized plasma modes, that smoothly interpolate between the relativistic and the non-relativistic regime. From the other side, we clarified the approximations implicit in the different approaches used in the small-momentum [14][15][16][17][18][19]34 or largemomentum 10-13,20,21 regime, showing that they can be obtained from the same gauge-invariant effective model once specific approximations are done. Finally, to make a direct connection with physically-meausurable quantities, we derived explicitly the density-density response of the layered superconductor, showing how the generalized plasmons are projected out in the physical response in the various energy/momentum regimes.
Our results, including the methodological aspects, are particularly relevant for future work aimed at investigating e.g. the plasma modes in confined geometries, as it is the case for surface plasmon polaritions 4,5 . The recent detection of surface Josephson plasma waves in an ultrathin film of the high-temperature superconductor La 1.85 Sr 0.15 CuO 4 shows that the field is now mature to apply near-field optics, efficiently used so far to investigate van der Waals layered metals 5 , to layered superconductors. At the same time, even if standard EELS or RIXS experiments have a rather low momentum resolution, with the lowest accessible momenta of about ∼ 0.1 of the reciprocal lattice unit, electron energy loss spectroscopy incorporated in a scanning transmission electron microscope (STEM-EELS) and equipped with a monochromator and aberration correctors has a high potential to combine high spatial and energy resolution 51,52 . As a consequence one can hope in the near future to be able to probe plasma excitations around the crossover scale k c ∼ 5µm −1 , where approximate solutions fails, and where both generalized plasma modes give a comparable contribution to the density response, as we showed in Sec. IV. In addition, at these wavevectors and energies one could explore the possible coupling with other collective modes like phonon, giving rise to hybrid phononpolaritons modes 4 , involving simultaneously transverse and longitudinal fields. As a final remark, it can be interesting to explore how the transverse-longitudinal mixing affects the plasmon dispersion in the normal state, even though in this case the soft Josephson plasmon is expected to be strongly overdamped by the quasiparticle continuum, especially in correlated systems like cuprates. Such a description can also help understanding how to model the Coulomb interaction among the scattered electrons and the local density of the layered metal in EELS experiments in cuprates, whose interpretation is still debated [46][47][48][49] . Indeed, a robust modelling of the plasma modes in conventional layered superconductors and metals is certainly a prerequisite in order to understand how charge fluctuations manifest instead in a correlated system like cuprates.
Φ ∆ (q) ≡ kσĉ + k−q/2σĉ k+q/2,σ being the pairing operator, U > 0 being the SC coupling constant and N denotes the number of lattice sites. In order to compute thermal averages over the hamiltonian (A1) we use the path integral formulation. Within such framework the the partition function of the system is given T and c, c are the Grassmann variables associated to the annihilation and creation operators respectively. To obtain the effective action in terms of the order-parameter collective degrees of freedom, the interacting action is decoupled in the particle-particle channel by means of the Hubbard-Stratonovich (HS) transformation by introducing the auxiliary complex field ∆: ∆(r) = (∆ 0 + δ∆(r))e θ(r) where ∆ 0 is the mean-field expectation value of the amplitude associated to the SC energy gap, δ∆ and θ are amplitude and phase fluctuations. By making an appropriate gauge transformation on the Grassmann fields c and c it is possible to make the dependence on the phase θ explicit in the action. One then finds that the HS transform of S I is phase-independent, while the free contribution S 0 becomes: where Ψ is the Nambu spinor, defined as the column vector c ↑ , c ↓ , and S BCS is the BCS saddle-point action, where only the mean-field value ∆ 0 is the complex field has been included, while fluctuations are contained inΣ. The BCS Green Function are defined from S BCS as: iω m = (2n + 1)πT being Matsubara fermionic frequencies, E k = ξ 2 k + ∆ 2 0 being the quasiparticles energy, τ i being the Pauli matrices.Σ is the self-energy, which depends, in principle, on both amplitude and phase fluctuations. Nonetheless, as long as one is interested in the low-temperature dynamics of phase fluctuations in layered cuprates, amplitude fluctuations can be neglected 12 . The self energy then reads: where ← → ∇ ≡ ∇ − ∇. Notice that, according to Goldstone Theorem, the phase θ appears only trough its time and spatial derivatives in Eq. (A7), i.e. no mass term for θ is allowed. Now, since the fermionic variables appear quadratically intoS 0 , one can integrate them out. Such procedure leads to the following effective action for phase fluctuations: where the trace is computed over both spin and momentum degrees of freedom. In order to study the phase dynamics we compute the effective action at Gaussian level, i.e. by including terms up to n = 2. It reads, in Fourier space: where q ≡ (iΩ m , q), Ω m = 2πmT being a bosonic Matsubara frequency, and: are the BCS response functions, which contain all the information on the microscopic fermionic degrees of freedom. Again, if one is interested in the low-temperature phase dynamics, the hydrodynamic limit of the gaussian action (A9) is the relevant one. Within such approximation, the BCS bubbles are computed in the static limit iΩ m = 0, q → 0, and one finds Eq. (10) of the main text. The formalism developed so far is appropriate for describing neutral superfluids. For the case of superconductors, the e.m. interaction between electrons must be taken into account. This can be achieved via the scalar and vector potentials φ and A, which account for internal e.m. degrees of freedom. They can be included into the self-energy (A7) by means of the minimal-coupling substitution, which reads in imaginary time as ∂ τ → ∂ τ −eφ, −i∇ → −i∇ + e c A for the fermionic degrees of freedom. After such a procedure the self energy (A7) undergoes the following modification: One immediately sees that Eq. (A13) can be obtained as well by making the substitution (7) on the phase degrees of freedom. At this point it is possible to expand the effective action at gaussian level in both the phase and the e.m. potentials: by doing such a procedure and by adding the free-e.m. contribution (11) one obtains the total action (13). External perturbations can be also introduced, in order to compute the response function of the system. For example, in order to compute the density-density response function, one can introduce a scalar perturbation δφ, which couples with the density operatorρ q = kσĉ + k+qσĉ kσ into the microscopic hamiltonian trough the following contribution: Eq. (A14) will give rise to the density insertion eδφτ 3 into the self energy (A13). At this point one can compute, trough the usual expansion (A8), the effective action at gaussian level in the phase, the internal e.m. degrees of freedom and the scalar perturbation. Then, once the g.i. variable ψ have been introduced and the phase has been integrated out, one is left with Eq. (50), which describes the coupling of the scalar perturbation δφ with the longitudinal degrees of freedom (described by ψ L ) of the system.