Collective modes and THz near field response of superconductors

We theoretically study the low energy electromagnetic response of BCS type superconductors focusing on propagating collective modes that are accessible with THz near field optics. The interesting frequency and momentum range is $\omega<\Delta$ and $q<1/\xi$ where $\Delta$ is the gap and $\xi$ is the coherence length. We show that it is possible to observe the superfluid plasmons, amplitude (Higgs) modes, Bardasis-Schrieffer modes and Carlson-Goldman modes using THz near field technique, although none of these modes couple linearly to far field radiation. Coupling of THz near field radiation to the amplitude mode requires particle-hole symmetry breaking while coupling to the Bardasis-Schrieffer mode does not and is typically stronger. The Carlson-Goldman mode appears in the near field reflection coefficient as a weak feature in the sub-THz frequency range. In a system of two superconducting layers with nanometer scale separation, an acoustic phase mode appears as the antisymmetric plasmon mode of the system. This mode leads to well defined resonance peaks in the near-field THz response and has strong anti crossings with the Bardasis-Schrieffer mode and amplitude mode, enhancing their response. In a slab of layered superconductor such as the high T$_c$ compounds, which can be viewed as a natural optical cavity, many branches of propagating Josephson plasmon modes couple to the THz near field radiation.


I. INTRODUCTION
Electromagnetic (EM) response is a fundamental property of superconductors. The response to static electric and magnetic fields (infinite conductivity and the Meissner effect) are the defining properties of the superconducting state. The response to the time dependent, very long wavelength transverse fields produced by far field radiation has been extensively studied [1][2][3][4][5] . Superconductors are also characterized by a diversity of sub-gap collective modes 6 including plasmons, acoustic phase modes, amplitude (Higgs) modes, the Carlson-Goldman modes, and the Bardasis-Schrieffer modes associated with fluctuations of subdominant order parameters, shown schematically in Fig. 1. For superconductors of current interest including cuprates 5,7 , iron pnictides 8 , NbSe 2 9 and MgB 2 10 the gap values and the relevant collective modes are in the terahertz (THz) range. These modes couple weakly, if at all, to far field transverse photons.
Recent progress in cryogenic near field nano optics 11 has enabled new generations of experiments probing the response of materials to short wavelength, primarily longitudinal, finite frequency electric fields [12][13][14] . This information is encoded in the near field reflection coefficient R p (ω, q) (see Appendix B). The essential new feature of the nano optics experiments as compared to traditional far field optics is the excitation of charge fluctuations in the material under study. In this paper we calculate the nano optics linear response R p (ω, q) of superconductors focusing on the contribution of collective modes. Each of the modes we consider couples to charge fluctuations and is therefore in principle observable in nano optics experiments. We calculate in detail the matrix elements coupling each mode to charge excitations and from this the signal of the nano optical response.
Charge fluctuations are constrained by the continuity equation, the proper treatment of which requires going be-yong BCS mean field theory in a manner consistent with the relevant Ward Identities [15][16][17][18] . In this paper we employ a one loop effective action method based on a Hubbard-Stratonovich transformation of the fundamental interacting electron system. This methodology is an efficient way to take order parameter and charge fluctuations into account while respecting conservation laws. It is equivalent to a diagrammatic calculation including vertex corrections 16 .
We now discuss the physics of the collective modes and FIG. 1. The schematic representation in the frequency-momentum plane of the collective modes that may appear in the electrodynamical response of a 2D superconductor. The blue area shows the low energy and long wave length region where weakly damped collective modes maybe observed. Anti crossing between the plasmon and Higgs mode and the BSM is not shown here. c is the speed of light, v F is the fermi velocity and D f is the normal state diffusion coefficient.
their coupling to light. The phase (Anderson-Bogoliubov-Goldstone) mode (heavy blue dashed line in Fig. 1) is the order parameter phase fluctuation. It is accompanied by a superfluid density oscillation 19 which in the presence of long range Coulomb interaction converts this mode into a plasmon 20 . In two dimensions (2D), the plasmon has a ω ∼ q dispersion, shown as the red line in Fig. 1 and discussed in section IV. The plasmons directly couple to near field radiation. In multi layer systems, mutual screening leads to branches of acoustic (linearly dispersing) plasmons (not shown in Fig. 1). These acoustic plasmons also couple to THz near field probes, see section VIII.
In the presence of abundant normal carriers as happens for example close to T c , the Coulomb potential of the superfluid density fluctuation can be screened and as a result, the gapless phase mode known as the Carlson-Goldman (CG) mode 21 (purple line in Fig. 1), appears. This mode is discussed in section VII.
The amplitude (Higgs) mode (green line in Fig. 1) is the gap amplitude fluctuation which couples to EM linear response only when particle hole symmetry is broken [22][23][24] and only at non-zero q because electron density fluctuation is needed to locally perturb the density of states and then the gap. The coupling is suppressed by the small parameter ∆/E F , as discussed in section V. Note that the Higgs mode does couple to far field radiation through third order nonlinear response which has been demonstrated experimentally [25][26][27] .
The Bardasis-Schrieffer mode (BSM) 28-31 is a fluctuation of subdominant pairing order parameters, e.g., d -wave fluctuations in a s-wave superconductor. Although it was proposed half a century ago 28 , the BSM has been very difficult to observe. Its signature was reported only recently in an iron based superconductor 32 . The BSM frequency is slightly below the gap for weak subdominant pairing and approaches zero as the subdominant pairing strength approaches the dominant one.
The rest of the paper is organized as follows. In section II we perform the Hubbard-Stratonovich transformation of the BCS Hamiltonian to obtain the Ginzburg-Landau effective action which includes the collective modes. Section III derives the linear EM response functions and their simple forms in the low energy limit. With the longitudinal optical conductivity, we analyse the properties of the collective modes in sections IV, V, VI and VII with a focus on 2D. Section VIII discusses the acoustic plasmon mode in superconducting double layer which is promising to be observed in THz near field optics. We then discuss the cluster of hyperbolic Josephson plasmons in naturally layered superconductors in section IX and show that they are greatly affected by the nonlocal correction to the optical conductivity and the discrete nature across the layers. Section X is a summary and conclusion, with pointers to the relevant equations and figures, for readers uninterested in the details of the derivations. Appendix A contains the definition and explicit forms of the correlation functions. Appendix B has the derivation of the reflection coefficients.

A. The action of fermion, gap and EM fields
The starting point is the BCS Lagrangian of attracting electrons coupled to electromagnetic (EM) field: is the electron annihilation operator, ξ(p) = ε(p) − µ and g < 0 is the attractive interaction. We will be interested in relatively low frequency longitudinal EM fluctuations where the magnetic field can be neglected. Thus d r 1 8π E 2 for 3D and → q 1 4π|q| E −q E q for a 2D plane embedded in 3D space where E q is the Fourier component of the electric field on the 2D plane. Note that the EM field A has an UV cutoff which is much smaller than the fermi momentum such that it mediates only the smooth part of the Coulomb interaction between the electrons. The high energy part of the photon is already integrated out which together with the phonons results in the effective interaction g whose cutoff is the Debye frequency ω D . For simplicity, we don't explicitly notate the photons in what follows.
Performing the Hubbard-Stratonovich transformation of the path integral Z = D[ψ, ψ]e −S in the pairing channel gives where the action describes coupled dynamics of the fermion field ψ, the EM field A and the gap ∆ l with l denoting the pairing angular momentum. Note that we have neglected detailed structure in g that is not important to our conclusions. The fermion kernel is where f l (p) describes the momentum dependence of the pairing function in each angular momentum channel. For simplicity, we consider only spin singlet pairing with the Nambu spinors being ψ † = (ψ † ↑ , ψ ↓ ). In the BCS regime of two dimensional superconductors, we can choose f l = cos(l θ k ) or sin(l θ k ) and the corresponding pairing interaction is g l = 1 2π d θ cos(θ)g (2k F sin(θ/2)). Note that for l = 0, the 1/2π factor should be changed to 1/4π. We assume that the l = 0 component of the interaction is the strongest and thus the ground state has s-wave pairing. We successively integrate out the fermion field to obtain the Ginzburg Landau action for the gap and EM field, and then the gap to obtain the action for EM field which gives the information of the EM response functions.

B. Integrating out the fermions
Expanding in the EM field, the Lagrangian density is where the Gor'kov Green's function for the Bogoliubov quasi particle is the paramagnetic contribution to the current is and the diamagnetic 'Drude' kernel is where the trace is over the frequency, momentum and spinor degrees of freedom of ψ.

C. Mean field as the saddle point
Assuming the ground state has s-wave pairing with ∆ independent on momentum, the mean field free energy is where E k = ξ 2 k + |∆| 2 = (ε k − µ) 2 + |∆| 2 is the quasi particle energy with gap ∆. Minimization of S mean field with respect to the gap renders the gap equation where f is the fermi distribution function. At zero temperature, the integral in Eq. (10) can be done and we obtain the condensation energy relative to the normal state: where ν is the density of state at the fermi level of the normal state. The gap at zero temperature is thus ∆ 0 = 2ω D e − 1 g ν for g ν 1. Without loss of generality, we pick the mean field gap ∆ to be real. The coherence length is defined as ξ = v F /∆. Note that the free energy is non analytic around ∆ = 0, i.e., the expansion coefficients in small ∆ all diverge. The Ginzburg Landau expansion in powers of ∆ is only possible at finite temperature and accurate close to T c .

D. Fluctuations
The fluctuations include those of the EM field, the s-wave gap and the subdominant pairing order parameters. The s-wave gap fluctuation can be separated into amplitude and phase: ∆ = (∆ 0 + ∆(r, t ))e i 2θ(r,t ) . It is convenient to perform a local gauge transformation 33 which changes the action Eq.
(3) to its equivalent form with the EM field changed to the gauge invariant ones In the effection action, the EM field always comes together with the phase fluctuations in the above gauge invariant form and therefore couples directly to the phase mode. The appearance of the other modes in the EM response can be inferred simply from their coupling to the phase mode. The action of fluctuations ∆ q around the mean field gap can be decomposed as is the phase action, is the amplitude action and is the action for the fluctuation of the sub dominant pairing channels, i.e., the Bardasis-Schrieffer modes. Finally, (19) is the coupling between phase and amplitude/BSM fluctuations. Global U (1) symmetry under θ → θ + δ is manifest here since the action depends only on derivatives of the phase. This ensures charge conservation.

E. Phase action
The quadratic kernel for the phase action is where χ (0) ρρ , χ (0) ρj and χ (0) jj are the density-density, densitycurrent and current-current correlation functions evaluated at the mean field saddle point. In the case of quadratic band, ε = p 2 /(2m), the system has Galilean invariance and D i j = n/mδ i j where n is the total carrier density.
Low energy limit-At zero temperature, for ω ∆ and ρj ∼ ωq and χ (0) jj ∼ q 2 , thus to leading order we have where ν is the normal state density of state at the fermi level. Therefore, the effective low energy Lagrangian of the phase fluctuation is 33 which describes the Nambu-Goldstone mode with velocity v g = n/(mν) = v F / d if EM field were not present, also known as the Anderson-Bogoliubov mode 2,15 . Here d is the space dimension. Due to the long range Coulomb interaction (coupling to EM field), the Goldstone mode does not actually exist but is shifted to the high frequency plasmons through the Anderson-Higgs mechanism 20 .

F. Amplitude action
The inverse propagator for amplitude fluctuation is At zero momentum q = 0 and rotated to real frequency, it has the analytical form where ∆ sc = ∆ and describes the quasi-particle effects. Specifically, F diverges as 1 2∆−ω as the frequency approaches 2∆ and has an imaginary part above 2∆ due to quasi-particle excitations. Thus G −1 a does not have a simple pole at ω = 2∆ and the Higgs amplitude mode is not well defined in the weak coupling BCS approximation 24 .
However, in systems with superconductivity coexisting with charge density wave (CDW), the quasi particle absorp- cd w is pushed beyond the Higgs frequency 2∆ SC and the latter becomes a well defined collective mode 23 . In this case, F is well behaved around the Higgs pole and we can approximate the propagator by where d is space dimension and the O(q 2 ) expansion can be found in Appendix A 1. Thus the Higgs mode frequency disperses roughly as ω 2 F q 2 as in Ref. 22 . Moreover, coupling between the Higgs mode and a higher frequency CDW phonon may further lower the Higgs mode frequency 22 as has been proposed for 2H -NbSe 2 9 .

G. BSM action
The inverse propagators for the fluctuations of the higher angular momentum pairing channels are where the correlator is defined in Appendix A 2. Note that there are two directions for the fluctuations of the subdominant order parameters on the complex plane: perpendicular to the mean field gap (in the 'imaginary' direction) and parallel to it (in the 'real' direction). The BSMs 28,30,31 are the 'imaginary' fluctuations, i.e., in the σ 2 channel. This channel of fluctuations have nonzero matrix element in exciting quasi particles on the gap edge, thus rendering χ σ 2 ,σ 2 (ω) diverging as the frequency approaches the gap 2∆ from below. As a result, the BSMs all have energies below 2∆. The 'real' modes, the fluctuations in the σ 1 channel, don't have poles and are not collective modes, as shown in Appendix A 2.
In this paper we focus on the d -wave BSM in an s-wave superconductor which couples to light already at quadratic order in q. In two dimension, there are two d -wave BSMs corresponding to d x 2 −y 2 and d x y . Assume the momentum is along x, to leading order in momentum, the inverse propagator of the d x 2 −y 2 mode is For g d ∈ (0, 2g ), at zero momentum, it has a pole below 2∆ which gives the mode frequency in the weak and strong BSM fluctuation limits. As g d changes from 0 to 2g , ω B S goes from 2∆ to 0. For g d > 2g , the ground state is no longer a s-wave one 29 . At finite momentum, the mode frequency shifts up as q 2 as shown in Appendix A 2.

H. Coupling terms
Under either time reversal or particle hole operation, the phase θ and BSM ∆ l are odd because they are fluctuations in the 'imaginary' direction on the complex plane while the amplitude fluctuation ∆ is even. Therefore, the linear coupling between phase and amplitude fluctuations breaks particle hole symmetry while that between phase and BSM does not.
Due to inversion and time reversal symmetry, the coupling coefficient can be expanded as where C 0 ,C i = 0 in a particle hole symmetric situation. Indeed these limits can be verified from their exact forms To study the linear EM response for q 1/ξ, it is sufficient to keep the leading terms. The simplest way to break the particle hole symmetry is to use an energy dependent electronic density of state (DOS), e.g., as in the parabolic band electron gas in three dimension. Assuming the DOS is g (ξ) = ν(1 + λξ/E F ) (note that ξ = ε − µ should not be confused with the coherence length), we obtain from and from χ (0) j∆ that see Appendix A 3 for detailed derivation. The small factor λ∆/E F characterizes the strength of particle hole symmetry broken.
In two dimension, from inversion, time reversal symmetry, that ∆ d being odd under time reversal and its d -wave character, the linear coupling coefficient between the phase and the d x 2 −y 2 BSM can be written as Diagrammatic representation of the EM linear response kernel Π, i.e., self energy of photon. First line is the photon self energy using the language of interacting electrons. The solid lines are electron Green's functions within the BCS approximation. The corresponding vertex correction should be included to restore the Ward identity 16 . Second line is the same thing but expressed in the language of coupling photons and order parameter fluctuations, as described by the Ginzburg-Landau action in Eq. (9). The first term is the bare current correlation K and the second term is the fluctuation contribution which corresponds to the vertex correction.
where B i = i π∆v 2 F F (ω) describes coupling between electric field and the BSM, see Appendix A 3. Since the angular momentum change is δl = 2 in exciting an s bound state to a d state, an inhomogeneous electric field is required to overcome the selection rule and thus the B i terms exist only at finite momentum. It will be shown that in the optical conductivity, at leading (O(q 2 )) order the B 0 term does not contribute and thus can be neglected. Note that coupling to BSM does not require breaking particle hole symmetry and is not suppressed by the small number ∆/E F . Thus in general, the BSM couples to EM more strongly than the Higgs mode.

A. Linear EM response
To obtain the linear EM response functions, one just needs to obtain the fluctuation action quadratic in the EM fields. Integrating out the amplitude results in with the kernel modified to For longitudinal EM response, it is inappropriate to directly employ the 'free' optical conductivity χ (0) jj +〈D i j 〉 or the density response χ (0) ρρ obtained from the BCS mean field Hamiltonian which breaks global U (1) gauge invariance. They actually do not satisfy charge conservation: K µν q ν = 0. The reason is that longitudinal EM fields excite order parameter phase fluctuations that are not captured by Bogoliubov quasi-particles.
The solution is to take into account the phase fluctuations which ensures charge conservation since the Euler-Lagrangian equation from Eq. (35) is just the continuity equation. By integrating out the phase (or simply solving the Euler Lagrangian equation of the phase), one finally obtains the EM action is the EM response tensor satisfying J µ = Π µν A ν and the continuity equation Π µν q ν = 0. Specifically, Π 00 = χ ρρ is the irreducible density density response (polarization function) and 1 i ω Π i j = σ i j is the optical conductivity. Note that 'irreducible' is with respect to Coulomb interaction, i.e., EM field. The diagrammatic representation of Eq. (38) is shown in Fig. 2. This is equivalent to correcting the current vertex by electron electron interactions after which gauge invariance 34 and thus Ward Identity 16 are recovered.
Eq. (38) should be viewed as the central result of this paper which contains all the information of linear coupling between EM field and the collective modes.
The major effect of disorder is to induce a nonzero σ 1 above twice the gap. The Mattis-Bardeen 1 theory for optical absorption completely relaxes momentum conservation in the quasi particle excitation process, and has proven accurate in various BCS type superconductors. In this paper, we implement the Mattis-Bardeen formula to describe the optical conductivity above the gap: where σ n is the normal state conductivity and E (x), K (x) are the complete elliptic integrals.

B. The low energy limit
At low temperature compared to T c , in the low energy limit of ω ∆ and q 1/ξ as shown in the blue region of Fig. 1, the electrodynamics can be described by the Lagrangian Eq. (22) which leads to the longitudinal optical conductivity: Here n s is the superfluid density and n s /m is the superfluid stiffness which in 3D is related to the magnetic penetration depth as λ B = c 2 4π m n s e 2 . This form closely resembles that of a hydrodynamic electron fluid 35,36 except that damping is completely suppressed here by the gap. Eq. (40) contains the whole crossover between the Drude limit ω v F q and the Thomas-Fermi limit ω v F q. Note that the amplitude/BSM and quasiparticle excitation don't enter here since they appear at higher energy. For a clean and isotropic BCS superconductor, v g = v F / d at zero temperature and gradually decreases to zero as temperature is raised to T c .
At non-zero temperature, one should add the contribution of the normal carriers which makes the conductivity into the 'two fluid' form For ∆ T , an analytical formula for the normal fluid conductivity with non-zero scattering rate can be found from the Boltzmann equation 37 . In the simple limits, where D n = πn n /m, γ is the scattering rate, is the normal state diffusion constant, d is space dimension and ν n is close to ν at temperate close to T c .
In principle, this two fluid formula Eq. (41) can be obtained from the general derivation Eq. (38) with electron-impurity or electron-phonon scattering taken into account. We sketch the derivation of the two fluid model here by calculating the polarization function from Eq. (38): Close to T c , we have where χ (0) come from the 'intra band' process among thermally excited quasiparticles while χ s is the 'interband' contribution from exciting quasi particle pairs which has the interpretation of superfluid susceptibility (compressibility). Close to T c , the χ (0) should resemble the polarization function of a normal fermi liquid, i.e., the Lindhard function 37,38 with non-zero scattering rate. Similarly 33,39 , where σ n is the 'intraband' part and is the superfluid density of a clean superconductor. Moreover, close to T c , the intraband contributions to the correlation functions should approximately satisfy the continuity equations ωχ (0) +qχ (0) ρj = 0 and q j σ ni j +ωχ (0) ρ j i = 0 since they are identical to those of the normal state at ∆ = 0. With these simplifications Eq. (43) becomes The corresponding conductivity is just Eq. (41) with σ n = χ (0) i ω/q 2 and v g = n s m /χ s .

IV. 2D PLASMONS
For simplicity, we neglect the coupling to the amplitude mode in this section. The plasmons are the charge density fluctuations and can be found by the zeros of the dielectric function where V q = 4π/q 2 for three dimension and V q = 2π/|q| for 2D. Together with Eq. (40), we obtain the plasmon dispersion ω p = 2D s q + v 2 g q 2 for two dimension. For three dimension, Eq. (48) predicts ω p = 4πn s e 2 /m + v 2 g q 2 ∆ which lies in the high energy regime beyond the limit of validity of our theory although the correct plasma frequency ω 2 p = 4πn s e 2 /m is obtained for a clean superconductor.
In the low frequency limit ω ω c , the 2D plasmon dispersion ω p = 2D s q approaches the edge of the continuum of vacuum propagating photons ω = c q (recall that here q is a two dimensional momentum and light modes disperse as ω = c q 2 + k 2 z for any k z ). For lower frequencies the analysis given here requires modification, because the electric fields associated with the plasmons begin to extend far from the 2D sheet, so that the plasmon couples much less strongly to near field radiation. The critical frequency can be estimated as ω c = 2D s /c corresponding to the energy ħω c = e 2 ħc ħ 2 n s m ∼ 1 137 E * F where E * F is the fermi energy equivalent to a two dimensional superfluid stiffness n s /m. For a clean, weakly correlated material E * F is of eV-scale and the crossover frequency is of the order of 1 THz. However, many superconductors of current interest 40 have much lower E * F so that the crossover frequency is well below the THz regime.
Assuming a doping level of n = 7×10 13 cm −2 (corresponding to the fermi momentum k F = 2π/(3 nm)) and the fermi velocity of v F = 2.5 × 10 5 m/s, one obtains a wavelength of λ ≈ 180 µm for the plasmon at 1 THz. This wavelength is close to that of the corresponding vacuum photon (300µm) although a substrate with large dielectric screening might make the plasmon wave length shorter. (a) At T = 0 K, the dominant feature comes from the plasmon while there is very weak anti crossing with the BSM at 3.0 THz. The coupling to Higgs mode is too weak to be seen. (b) At T = 79.6 K close to T c = 80 K, the plasmon is overdamped while the CG mode appears as a weak crossover of R p . Note the difference in color scales between (a) and (b). The fermi momentum is k F = 2π/(3 nm), the fermi velocity is v F = 2.5 × 10 5 m/s, the normal state scattering rate is γ = 30 THz, the gap at zero temperature is ∆ = 3.0 THz and Nevertheless, in a dirty superconductor with a large normal state scattering rate γ T c , the sub-gap plasmon frequency is mainly determined by the superfluid density which is only a part of the total density even at zero temperature: n s ∼ nT c /γ. At the same THz frequency far below γ, the plasmon wavelength is smaller by the factor T c /γ and they become more confined to the 2D plane. Below the gap, weakly damped plasmons with dispersion ω = 2D s q couple strongly to near field probe as shown by the near field reflection coefficient in Fig. 3 (a). For the ratio T c /γ = 0.1, the 1 THz plasmon wave length is shrunk by the same factor to 18 µm.

V. HIGGS MODE
The Higgs mode couples to phase fluctuation in Eq. (19), manifests itself in the second term of Eq. (36) and finally enters the EM response through Eq. (38). The density response with the Higgs mode correction is thus Since the C 0 and C i terms contribute terms at the same order, we take C i = 0 to arrive at a simplified expression Thus the longitudinal optical conductivity is where is the dimensionless coupling constant of the Higgs mode to EM and ω hq = 4∆ 2 sc + v 2 F q 2 /d is the Higgs mode frequency. Since λ is order one and sinh −1 ω D ∆ is not a large number, this coupling is simply suppressed by the small number ∆ E F . Note that at q = 0, the optical conductivity reduces to the Drude form and there is no signature of Higgs mode in it. Therefore it is impossible to observe it in far field THz linear response. In contrast, near field optical imaging technique has access to non-zero q where the Higgs mode manifests itself through coupling to the plasmons.
Specifically, for a monolayer superconductor, the coupled collective modes can be found as the poles of Eq. (48). The weight of the Higgs pole in R p scales as W hi g g s ∼ κ 2 v 2 g q 2 /∆ for ω h ω p , i.e., well before the Higgs mode crosses the plasmon. Nevertheless, the most prominent signature of the Higgs mode is its anti crossing with the plasmon mode which happens roughly at ω p (q) = ω hq . A detailed solution of Eq. (48) gives the frequency splitting at the anti-crossing as where q is the momentum at the anti-crossing. Therefore, the splitting will be bigger if the anti-crossing happens at larger momentum.

VI. BARDASIS-SCHRIEFFER MODE
Optical excitation of the BSM can be viewed as transition from an s bound state of the cooper pair to a d bound state. This is forbidden in far field optics for two reasons: first, unlike the Hydrogen atom case, uniform electric field exerts the same force on the two electrons and does not change the internal structure; second, both s-state and d-state have even parity which forbids the transition due to optical selection rule. Thus it is necessary to go to nonzero momentum for its nonzero coupling to EM field. Indeed, the coupling constant is proportional to ξq which is appreciable when the electric field becomes substantially nonuniform on the scale of a cooper pair size.
Plugging Eq. (34) into Eq. (36) and Eq. (38) gives the appearance of the BSM in the longitudinal optical conductivity: where κ B S = π 8 v 2 F /v 2 g ∼ 1 is the dimensionless coupling constant between BSM and EM. The B 0 terms are higher order in q and are neglected. Note that if the momentum q is along x, the BSM means the d x 2 −y 2 order parameter fluctuation.
The BSM couples to near field more strongly than the Higgs mode due to the absence of the ∆/E F factor in the coupling constant κ B S . Solving the pole equation, Eq. (48), one obtains the frequency splitting at the anti-crossing between BS and plasmon δω ≈ κ B S v g q (56) which scales linearly with the momentum at the anticrossing.

VII. CARLSON-GOLDMAN MODE
The CG mode is a superfluid density fluctuation accompanied by the counter flow of normal carriers such that the Coulomb potential from the superfluid fluctuation is almost completely screened 21,[41][42][43][44] . This screening requires a large density of normal carriers which is typically found near T c . The speed v g of the CG mode depends on the ratio between the superfluid density n s and superfluid susceptibility χ s = π 4 ∆ T c ν and has different expressions in the clean and dirty limits: where ζ(x) is the Riemann Zeta function and we have used the fact that n s = n π∆ 2 2γT c for dirty superconductors and n s = 2(1 − T /T c )n for clean superconductors close to T c .
Its dispersion can be derived from the two fluid conductivity Eq. (41) by setting = 0 which yields Note that the CG mode can be understood as a sound with the standard sound velocity n s m /χ s and χ s being the superfluid compressibility. The latter is smaller than ν, the compressibility of the whole fluid in the low frequency thermal dynamic limit, because the super and normal fluids move out of phase in this relative high frequency regime. The local accumulation of superfluid causes the local chemical potential to shift up leading to change of quasiparticle energy and charge. However, the quasiparticle occupation number relaxes too slowly and cannot adjust itself to this change 42 , resulting in 'branch imbalance' 46 as shown in FIG. 5. Physical picture of the local chemical potential shift and the quasi particle occupation in the CG mode. The plus and minus signs indicate the signs of the quasi particle charge. This quasiparticle distribution is referred to as 'branch imbalance' or 'charge mode' in the literature 46 . Fig. 5. Normal impurities cannot relax this branch imbalance because in s-wave superconductors, the elastic scattering matrix element u k u k − v k v k vanishes between hole like and electron like states at the same energy. Inelastic scattering due to, e.g., phonons, does relax branch imbalance and cause extra damping to the CG mode but we assume it to be small. In d -wave superconductors, the same matrix element is non-zero due to anisotropy of the gap which allows normal impurities to relax the branch imbalance and bring extra damping to the CG mode 47 .
In clean superconductors the CG mode can cross the diffusion line before reaching the gap, entering the regime ω D f q 2 where the normal fluid part of the Eq. (41) is in the Thomas-Fermi form Eq. (42). The normal fluid still screens the CG mode but with a Thomas-Fermi screening character. The CG mode speed is slightly modified to v CG = v 2 g + n s ν n m in this regime but still remains close to v g since the second term is much smaller.
The original experiment using Josephson tunneling junctions by Carlson and Goldman 21 seems to be the only observation of this novel collective mode, the closest analogy to the Goldstone mode of U (1) symmetry breaking in a superconductor. In the optical conductivity measured by far field optics, the CG mode might move part of the superfluid spectra weight to finite frequency due to smooth disorder 48,49 . At non-zero momentum, being almost charge neutral, the CG mode appears as a very weak feature in the near field reflection coefficient: a one percent crossover of Abs[R p ] as shown by Fig. 3(b) plotted for a typical dirty superconductor close to its T c .

VIII. DOUBLE LAYER SUPERCONDUCTOR
Consider the system made of two superconducting layers separated by some small distance a, as shown in Fig. 6. We neglect the Josephson coupling between the layers which is weak for a substantially larger than atomic scale. Each layer has an in plane conductivity described by Eq. (40) at low temperature. In the quasi static limit, the 2D plasmon dispersion can be obtained from the following eigenmode condition which leads to two plasmon branches The upper branch is the symmetric mode whose dispersion follow the ω + ∼ q law at small momentum. The lower (anti symmetric) branch is an acoustic mode which has the dispersion for q 1/a. This acoustic mode is charge fluctuations of the two layers which are out of phase such that the net charge fluctuation is near zero if looked at far away. In other words, the Coulomb interaction is mutually screened and is modified to the effective short range form V (q) = 2π(1 − e −aq )/q that makes the mode acoustic. Both modes correspond to non-zero momentum oscillations of the phase of the superconducting order parameter. This acoustic plasmon can be viewed as the Goldstone mode which recovers its acoustic nature because Coulomb interaction is greatly weakened. Its speed still has a large contribution 2D a ∼ αk F av F from the residual Coulomb interaction where α = e 2 /(ħv F ) is the 'fine structure constant'. In BSCCO 2212 at typical doping 50 , α ≈ 9 since v F ≈ 2.5 × 10 5 m/s. For k F = 2π/(10 nm) and a = 3 nm, the ratio between the speeds of this acoustic mode and the original Goldstone mode is v − /v g ≈ 6 which means they are at the same order of magnitude. Therefore, an accurate measurement of the acoustic plasmon dispersion would contain the information of the 'Goldstone mode' speed.
In order for the acoustic mode to be observable to near field probe, it should have substantial spectral weight in the the near field reflection coefficient derived in Appendix B. Given the same amplitude of charge density oscillation in each layer, the electric field generated by the two layers tend to cancel each other since they are opposite in sign. The remaining field is weaker than the symmetric plasmon mode by a factor of q a/2 and the near field spectra weight is weaker by (q a/2) 3/2 . Nevertheless, the acoustic mode is still visible as shown by the R p plotted in Fig. 6. Moreover, since the acoustic plasmon has higher momentum given the same frequency in the THz range, it has stronger coupling to the Higgs/BSMs. Thus there is more prominent anti crossing feature between them, as shown in Fig. 6. Note that Eq. (54) and Eq. (56) apply to anti-crossings with both the symmetric and anti-symmetric modes. For example, the anti crossing of the BSM with the acoustic plasmon happens at a momentum roughly 20 times that with the symmetric plasmon, rendering the energy splitting 20 times larger than the latter.

IX. BULK LAYERED SUPERCONDUCTORS
For the class of layered superconductors, e.g., High T c cuprates, there is Josephson coupling between the layers and the low temperature and subgap collective modes are the Josephson plasmons 5 . Considering only the phase degree of freedom, the Lagrangian for an evenly spacing layered superconductor is Ad z where θ n (r ), φ n (r ) and A n (r ) are the phase, scalar and vector potentials on the nth layer and E c is the Josephson coupling energy per unit area and we have set e = 1. For longitudinal fields we are interested in, we can choose tha gauge where A n (r ) = 0. Due to continuous translational symmetry in plane and discrete in z direction, it is convenient to Fourier transform the fields into the 'Bloch' form where a is the layer spacing, q is the in plane momentum and k z ∈ (−π/a, π/a) is the lattice momentum in z direction. The Lagrangian Eq. (65) is diagonalized as Solving the Euler-Lagrange equation of the phase and making use of the expression of the charge density ρ = ν(∂ t θ−φ), we obtain the 'nonlocal' polarization function Due to the discrete Fourier transform, the Coulomb potential kernel is modified to V (k z , q) = 2πe 2 q sinh(aq) cosh(aq) − cos(ak z ) . the Carlson-Goldman mode appears but as a very weak feature across the resonance since it has almost no net charge density fluctuation. Note that this mode is not enhanced in multilayer systems. The amplitude (Higgs) mode appears in the EM response because it couples to the phase fluctuation with a matrix element that is non-zero if there is no perfect particle hole symmetry (Eq. (32)). The ultimate coupling to THz near field is proportional to the square of the near field momentum q (Eq. (52)), and is strongly enhanced by an anticrossing with the plasmon or phase modes. The Higgs mode is only weakly visible for monolayer materials because the q plasmon dispersion means that the anticrossing occurs at a very small momentum (Fig. 3(a)). The feature is much more easily visible in bilayer systems as an anti crossing with the acoustic plasmon (phase) mode ( Fig. 6(b)). The coupling to Bardasis-Schrieffer (subdominant order parameter) mode is very similar to that of the Higgs mode, except that it does not require particle hole symmetry breaking. It is again most easily visible as a large q anti-crossing with the phase (Fig. 6(b)) or plasmon mode (Fig. 3(a)).
In multilayer superconductors, a multiplicity of phase modes exist, coined the hyperbolic Josephson plasmons (Fig. 7). The plasmon dispersion is hyperbolic ( < 0 for in plane and > 0 for out of plane), leading to total internal reflection ( Fig. 7(a)) and many plasmon branches with Higgs and BSMs visible as anti-crossings. The multiplayer nature means there are multiple branches of Higgs modes/BSMs, but they are weakly separated and may be difficult to resolve.
On the experimental side, detection of the collective modes offers useful information about both the ground state and the low lying excited states. On the theory side, knowledge of how to excite the collective modes are often the first step towards understanding non equilibrium dynamics. From the technological point of view, the low loss plasmonic modes are promising as information carriers in superconductor wave guides. The multiplayer systems described in Sections VIII and IX can be viewed as a kind of naturally occuring photonic cavities which enhance light matter coupling.
The formalism presented here is for s-wave superconductors. For d -wave superconductors, the qualitative features of the EM response such as the two fluid model Eq. (41), the plasmons and the Carlson-Goldman mode should be the same. Nevertheless, the CG mode might exist down to much lower temperature because of the large proportion of the normal fluid 56 although it might be heavily damped by normal disorder 47 . Due to the nodes in the d -wave gap, the THz plasmons might experience substantial damping even at zero temperature. The effect of disorder is not explicitly taken into account and would be a useful extension of the present research, e.g., disorder assisted Cherenkov radiation of plasmons by quasiparticles. is the electron Green's function. Rotation from imaginary to real time makes the time ordered correlation functions into retarded ones. In the correlation functions involving the currents, one should change the σ vertex to the current vertex. For example, Evaluating the correlation function Eq. (A1) renders With the knowledge of the gap equation, the Higgs propagator is thus 24 where is shown in Fig. 9. At non-zero momentum, the propagator is 22 whose O(q 2 ) expansion gives (A11)

The BS propagator
The total order parameter can be written as ∆ k = ∆ + l ∆ l (r, t ) f l (k) where we have chosen the mean field gap ∆ to be real. The subdominant pairing order parameter fluctuations ∆ l can have two possible directions: 1, orthogonal to ∆ on the complex plane or in the 'imaginary' direction; 2, parallel to ∆ or in the 'real' direction. The 'imaginary' fluctuations are the BSMs while the 'real' ones don't have poles and are not collective modes.

The linear coupling between phase and BSM/Higgs modes
The coupling between phase fluctuation and the Higgs mode requires particle hole symmetry breaking which we model using an energy dependent DOS g (ξ) = ν(1 + λξ/E F ). The coupling constants are derived from the correlation functions in Eq. (31). From the general formula Eq. (A6) at zero temperature, the temporal part is Written in matrix form, Eq. (B2) becomes (B3) whose solution gives the linear relation between the fields each side of the 2D layer whereM is the transfer matrix. Setting φ 2↑ = 0, one obtains the near field reflection coefficient for a 2D layer ω σ is the dielectric function in 2D.

Double layer
As shown in Fig. 13(b), applying the reflection problem twice, one obtains Setting φ 3↑ = 0 yields the reflection coefficient Eq. (64) for the double layer system. A characteristic plot of the reflection coefficient of the double layer system is Fig. 12 where resonances due to the symmetric and anti symmetric plasmons show up.

A slab with nonlocal optical response
If the polarization function is nonlocal, more unfortunately, if it depends also on the z direction momentum k z such as that of the layered superconductor Eq. (68), the near field reflection coefficient of the vacuum-infinite superconductor interface should be modified to R p (ω, q) = i q − k z z (ω, q, k z ) i q + k z z (ω, q, k z ) where q is the in-plane momentum which is a conserved quantity and k z is that in the nonlocal medium determined by the condition (ω, q, k z ) = 0. For a slab with finite thickness a, as shown in Fig. 13(c), the transfer matrix method for solving the reflection problem renders R sl ab (ω, q) = R p 1 − e 2i k z d where R p is from Eq. (B7).