{'fraction_non_alphanumeric': 0.046763206375048195, 'fraction_numerical': 0.017222912471616468, 'mean_word_length': 4.044084278768233, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enhancement of the electron-phonon scattering induced by intrinsic surface plasmon-phonon polaritons

We investigate light-matter coupling in metallic crystals where plasmons coexist with phonons exhibiting large oscillator strength. We demonstrate theoretically that this coexistence can lead to strong light-matter interactions without external resonators. When the frequencies of plasmons and phonons are comparable, hybridization of these collective matter modes occurs in the crystal. We show that the coupling of these modes to photonic degrees of freedom gives rise to intrinsic surface plasmon-phonon polaritons, which offer the unique possibility to control the phonon properties by tuning the electron density and the crystal thickness. In particular, dressed phonons with reduced frequency and large wave vectors arise in the case of quasi-2D crystals, which leads to large enhancements of the electron-phonon scattering in the vibrational ultrastrong coupling regime. This suggests that photons can play a key role in determining the quantum properties of certain materials. A non-perturbative self-consistent Hamiltonian method is presented that is valid for arbitrarily large coupling strengths.

Using light to shape the properties of quantum materials is a long-standing goal in physics, which is still attracting much attention. [1][2][3] Since the 60s, it is known that superconductivity can be stimulated via an amplification of the gap induced by an external electromagnetic radiation. 4,5 Ultrafast pump-probe techniques have been recently used to control the phases of materials such as magnetoresistive manganites, 6 layered high-temperature superconductors, 7 as well as alkali-doped fullerenes, 8 by tuning to a specific mid-infrared molecular vibration. It is an interesting question whether the phases of quantum materials may be also passively modified by strong light-matter interactions in the steady-state, without external radiation. Recently, the possibility of modifying superconductivity by coupling a vibrational mode to the vacuum field of a cavity-type structure has been suggested. 9 This idea was then investigated theoretically in the case of a FeSe/SrTiO 3 superconducting heterostructure embedded in a Fabry-Perot cavity. 10 In this work, enhanced superconductivity relies, however, on unrealistically large values of the phonon-photon coupling strength, in a regime where light and matter degrees of freedom totally decouple. 11 This occurs beyond the ultrastrong coupling regime, which is defined when the light-matter coupling strength reaches a few tens of percents of the relevant transition frequency. [12][13][14] Thanks to the confinement of light below the diffraction limit, 15 an alternative approach to engineer strong [16][17][18][19][20][21][22][23][24][25][26] and even ultrastrong light-matter coupling [27][28][29] is the use of plasmonic resonators. The latter allow the propagation of surface plasmon polaritons (SPPs), which are evanescent waves originating from the coupling between light and collective electronic modes in metals, called plasmons. 30 Interestingly, most quantized models involving SPPs are well suited to describe the strong coupling regime, 24,31 while a corresponding theory for the ultrastrong coupling regime is still missing.
In contrast to SPPs, an Hamiltonian model valid for all interaction regimes of surface phonon polaritons stemming from the coupling between light and ionic charges already exists. 32 In the absence of external resonator, strong coupling between SPPs in a 2D material and surface phonons of the substrate has recently stimulated considerable interest. [33][34][35][36][37] Such a situation can not generally occur in the same crystal due to the screening of the ion charges by electrons. Nevertheless, exceptions exist as a result of the coupling between phonons and interband transitions. 38,39 These include bilayer graphene, 40 alkali fullerides A 6 C 60 , 41 organic conductors (e.g. K − TCNQ 42,43 ), as well as some transition metal compounds. 44 Interestingly, some of these materials or similar compounds feature superconductivity and charge density waves under certain conditions, which are believed to be driven, at least partly, by electron-phonon interactions. [45][46][47][48][49] It is an open question whether the ultrastrong coupling stemming from the coexistence of plasmons and phonons with large oscillator strengths can affect the electron-phonon scattering without external resonator.
In this work, we propose and investigate the possibility to tune material properties of certain crystals via an hybridization of phonons, plasmons, and photons which is intrinsic to the material, and without the use of an external resonator. The coexistence of plasmons and phonons within the same crystal gives rise to hybrid plasmon-phonon modes, which have been extensively studied in weakly-doped semiconductors and semimetals. [50][51][52][53][54] Here, we show that while hybrid plasmon-phonon modes do not affect the electron-phonon scattering at the level of the random phase approximation, strong interactions between these modes and photons offer a unique possibility to control the energy and momentum of the resulting surface plasmon-phonon polaritons by tuning the electron density and the crystal thickness: "Dressed" phonons with reduced frequency arise in the ultrastrong coupling regime, where the electronic and ionic plasma frequencies become comparable to the phonon frequency.
These dressed phonons are shown to exhibit unusually large momenta comparable to the Fermi wave vector in the case of quasi-2D crystals, which leads to large enhancements of the electron-phonon scattering. These results suggest that photons can play a key role in determining the quantum properties of certain materials.
We utilize a self-consistent Hamiltonian method based on a generalization of that introduced in Ref. [55], which is valid for all regimes of interactions including the ultrastrong coupling regime, and in the absence of dissipation. This quantum description of surface plasmon-phonon polaritons provides an ideal framework to investigate how the latter can affect the quantum properties of the crystal. Furthermore, our method can be generalized to various surrounding media and geometries of interest such as layered metallo-dielectric metamaterials and simple nano-structures, and differs from the usual effective quantum description of strong coupling between excitons and quantized SPPs. 24,31 We show in the supplemental material that the latter model leads to unphysical behaviors in the ultrastrong coupling regime.
The paper is organized as follows: i) We first introduce the total Hamiltonian of the system under consideration; and ii) Diagonalize the matter part leading to plasmon-phonon hybridization in the crystal; iii) We derive the coupling Hamiltonian of these hybrid modes to photonic degrees of freedom, which gives rise to intrinsic surface plasmon-phonon polaritons. iv) We determine the Hamiltonian parameters and its eigenvalues self-consistently using the Helmoltz equation, and characterize the properties of the resulting surface plasmon-phonon polaritons; v) We show that while plasmon-phonon hybridization can not solely modify the electron-phonon scattering at the level of the random phase approximation (RPA), the coupling of these hybrid modes to photons can lead to large enhancements of the electron-phonon scattering in the crystal.

Quantum Hamiltonian
We consider a metallic crystal of surface S and thickness in air, which contains a free electron gas and a transverse optical phonon mode with frequency ω 0 . The excitation spectrum of the electron gas features a collective, long-wavelength plasmon mode with plasma frequency ω pl , as well as a continuum of individual excitations called electronic dark modes that are orthogonal to the plasmon mode. Both plasmons and phonons are polarized in the directions u z and u // (Fig. 1). The Hamiltonian is derived in the Power-Zienau-Woolley representation, 56 which ensures proper inclusion of all photon-mediated dipole-dipole inter-actions, and can be decomposed as H = H pol + H el−pn with H pol = H pt + H mat−pt + H mat .
The first term in H pol reads H pt = 1 2 0 dR D 2 (R) + 1 2 0 c 2 dR H 2 (R), and corresponds to the photon Hamiltonian. c denotes the speed of light in vacuum, R ≡ (r, z) the 3D position, 0 the vacuum permittivity, and D and H are the displacement and magnetic fields, respectively. The light-matter coupling term reads H mat−pt = − 1 0 dR P(R) · D(R). Here, P = P pl +P pn denotes the matter polarization field associated to the dipole moment density, where P pl and P pn correspond respectively to the plasmon and phonon contributions 1 . The matter Hamiltonian is decomposed as H mat = H pn + H pl + H P 2 , where H pn and H pl denote the contributions of free phonons and plasmons which are provided in the supplemental material, and H P 2 = 1 2 0 dR P 2 (R) contains terms ∝ P 2 pl , P 2 ph , and a direct plasmon-phonon interaction ∝ P pl · P ph . Finally, H el−pn includes the contribution of the free, individual electrons (electronic dark modes), as well as the coupling Hamiltonian between electronic dark modes and phonons.
Air Air 0 0.5 1 1.5 2 0 1 2 Figure 1: Intrinsic surface plasmon-phonon polaritons: Metal of thickness in air featuring a plasmon mode with frequency ω pl (gray dashed line) and an optical phonon mode with frequency ω 0 (green dashed line). Both plasmons and phonons are polarized in the directions u z and u // . When ν pn = 0 and ω pl ∼ ω 0 , the matter Hamiltonian H mat provides two hybrid plasmon-phonon modes with frequencies ω 1 and ω 2 [Eq. (1)] represented as a function of ω pl /ω 0 for ν pn = 0.5ω 0 . The surface polaritons generated by the coupling of these modes to light propagate along the two metal-dielectric interfaces with in-plane wave vector q.

Diagonalization of the matter Hamiltonian
Hybrid plasmon-phonon modes have been extensively studied in weakly-doped semiconductors and semimetals. [50][51][52][53][54] In this section, we propose a derivation of these hybrid modes using a non-perturbative quantum method. In order to diagonalize H mat , the phonon polarization P pn is written in terms of the bosonic phonon annihilation and creation operators B Q,α and B † Q,α , with the 3D wave vector Q ≡ (q, q z ) and α = //, z the polarization index. The lattice polarization field P pn is proportional to the ion plasma frequency ν pn = Z 2 N M 0 a 3 , which plays the role of the phonon-photon coupling strength. Here, N denotes the number of vibrating ions with effective mass M and charge Z in a unit cell of volume a 3 . The plasmon polarization P pl is provided by the dipolar description of the free electron gas corresponding to that of the RPA. 55 In the long-wavelength regime, P pl is written in terms of the bosonic plasmon operators P KQ and P † KQ (K denotes the 3D electron wave vector), superpositions of electron-hole excitations with wave vector Q across the Fermi surface. As explained in the supplemental material, H mat can be put in the diagonal form H mat = Q,α,j=1,2 ω j Π † Qjα Π Qjα , where the hybrid plasmon-phonon operators Π Qjα are superpositions of B Q,α , P KQ , and their hermitian conjugates, and satisfy the bosonic commutation relations [Π Qjα , Π † Q j α ] = δ Q,Q δ j,j δ α,α . The hybrid mode frequencies ω j are represented on Fig. 1 and are given by Here, j = 1 and j = 2 refer respectively to the signs − and +, and the longitudinal phonon mode with frequency ω 0 = ω 2 0 + ν 2 pn is determined by the combination of the short-range restoring force related to the transverse phonon resonance, and the long-range Coulomb force associated to the ion plasma frequency. The diagonalization procedure of H mat further provides the well-known transverse dielectric function of the crystal: 50 One can then use the eigenmodes basis of the matter Hamiltonian to express the total polarization field as with frequency ω 0 in the case ω pl ω 0 (Fig. 1). In the following, we focus on the most interesting case occuring close to the resonance ω pl ∼ ω 0 , where plasmon-phonon hybridization occurs.

Coupling to photons
Due to the breaking of translational invariance in the z direction, the 3D photon wave vector can be split into an in-plane and a transverse component in each media as Q = qu // + iγ n u z , where n = d, cr refers to the dielectric medium and the crystal, respectively.
For a given interface lying at the height z 0 , the electromagnetic field associated to surface waves decays exponentially on both sides as e ±γn(z−z 0 ) with the (real) penetration depth γ n . The displacement and magnetic fields can be written as superpositions of the fields generated by each interfaces m = 1, 2, namely D(R) = q,m 4 0 c S e iq·r u qm (z)D qm and H(R) = q,m w q 4 0 c S e iq·r v qm (z)H qm , with w q the frequency of the surface waves (still undertermined) and q ≡ |q|. The mode profile functions u qm (z) and v qm (z) depending on γ n are provided in the supplemental material. Using these expressions, the photon Hamiltonian takes the form where the overlap matrix elements A mm q and B mm q depend on the parameters γ cr and γ d .
The next step corresponds to finding the electromagnetic field eigenmodes, which consist of a symmetric and an antisymmetric mode, such that the photon Hamiltonian can be put in the diagonal form: C qσ is determined using the Ampère's circuital law, which provides C qσ = wq 2cβqσ . 55 The light-matter coupling Hamiltonian H mat−pt is derived using the expression of D(R) in the new basis together with Eq. (3). While q is a good quantum number due to the in-plane translational invariance, the perpendicular wave vector of the 3D matter modes q z is not. In the light-matter coupling Hamiltonian H mat−pt , photon modes with a given q interact with linear superpositions of the 3D matter modes exhibiting different q z . The latter are denoted as quasi-2D "bright" modes, and are defined as π qσj = qz,α f ασ (Q)Π Qjα , where f ασ (Q) stem from the overlap between the displacement and the polarization fields and is determined by imposing the commutation relations [π qσj , π † q σ j ] = δ q,q δ σ,σ δ j,j .
Using a unitary transformation (see supplemental material), the matter Hamiltonian can be decomposed as H mat = q,σ,j ω j π † qσj π qσj + H dark , where the second term is the contribution of the "dark" modes that are orthogonal to the bright ones, and which do not interact with photons. Without this contribution, the polariton Hamiltonian reads H pol = q,σ H qσ with and the vacuum Rabi frequency

Surface plasmon-phonon polaritons
The Hamiltonian H pol exhibits three eigenvalues in each subspace (q, σ), and only the lowest two (refered to as lower and upper polaritons) correspond to surface modes (below the light cone). At this point, we have built an Hamiltonian theory providing a relation between the field penetration depths γ n and the surface wave frequencies w qσζ , where ζ = LP, UP refers to the lower (LP) and upper (UP) polaritons. This eigenvalue equation can be combined with the Helmholtz equation n (ω)ω 2 /c 2 = q 2 − γ 2 n in order to determine these parameters. 55 We use a self-consistent algorithm which starts with a given frequency w qσLP , then determine γ n from the Helmholtz equation with d = 1 and cr (w qσLP ) given by Eq.
(2), and use these γ n to compute the parameters entering the Hamiltonian H pol . The latter is diagonalized numerically (see supplemental material), which allows to determine the new w qσLP . The algorithm is applied independently for the symmetric (σ = +) and antisymmetric (σ = −) modes until convergence, which is ensured by the discontinuity of both light and matter fields at each interface. 55 We now use this method to study the surface polaritons in our system. As an example, we consider the case ω pl = 1.5ω 0 with different crystal thickness q 0 = 10 [ Fig. 2 a)] and q 0 = 0.01 [ Fig. 2 b)], and compute the surface polariton frequencies w qσζ (top panels), as well as their phonon (i = pn), plasmon (i = pl), and photon (i = pt) admixtures W LP i,qσ for ν pn = 0 and ν pn = 0.5ω 0 (bottom panels). Precise definitions of these quantities are provided in the supplemental material, and and q are both normalized to q 0 = ω 0 /c. Considering typical mid-infrared phonons with ω 0 ∼ 0.2eV, the two dimensionless parameters q 0 = 10 and q 0 = 0.01 correspond to = 10µm (first case) and ∼ 10nm (second case), respectively.
In the first case [ Fig. 2 a)], the two surface modes at each interface have negligible overlap and the modes σ ± therefore coincide. For ν pn = 0, the plasmon-photon coupling is responsible for the appearance of a SPP mode (black line) with frequency w 0 q , which enters in resonance with the phonon mode at q r ≈ 2.2q 0 . While for q q 0 this SPP is mainly composed of light (w 0 q ∼ qc), it features an hybrid plasmon-photon character for q ∼ q 0 , and becomes mostly plasmon-like as w 0 q approaches the surface plasmon frequency ω pl / √ 2 for q q 0 . For ν pn = 0, a splitting between the two polaritons branches (thick red and blue lines) is clearly visible, and the latter consist of a mix between phonons, plasmons, and photons in the vicinity of q = q 0 . In the regime q > q 0 , since the mostly phonon-like LP exhibits a ∼ 25% plasmon admixture, this polariton mode can be seen as a "dressed" phonon with frequency red-shifted from ω 0 . Similarly, the "dressed" plasmon mode (UP) is blue-shifted from the surface mode frequency ω pl / √ 2 due its ∼ 25% phonon weight.
The frequencies of the LP (red line) and UP (blue line) at resonance are represented as a function of ν pn /ω 0 on the inset. Here, the resonance is defined by the condition w 0 q = ω 0 , which provides q r ≈ 2.2q 0 . We observe that the polariton splitting w qσUP −w qσLP is symmetric with respect to w qσζ = ω 0 (black dotted line) for ν pn ω 0 , and becomes slightly asymmetric in the ultrastrong coupling regime, 12 when ν pn is a non-negligible fraction of ω 0 . Furthermore, the splitting is found to decrease rapidly as the electronic plasma frequency ω pl is increased.
In the second case q 0 = 0.01 [ Fig. 2 b)], the crystal is thinner than the penetration depth γ cr of the surface waves at the two interfaces, and the latter overlap. This results in two sets of modes σ = ± with different frequencies. For ν pn = 0, the symmetric and antisymmetric SPPs with frequency w 0 qσ=± are represented as black dashed and solid lines, respectively.
The resonance between the symmetric SPP and the phonon mode occurs for q ∼ q Resonance between SPPs and phonons (horizontal dotted lines) is indicated by thin vertical lines. The antisymmetric (σ = −) and symmetric (σ = +) modes correspond to the solid and dashed lines, and the light cone with boundary ω = qc is represented a grey shaded region. Inset: Frequency w qσζ /ω 0 of the lower (red line) and upper (blue line) surface polaritons versus ν pn /ω 0 for q 0 = 10 and ω pl = 1.5ω 0 . The phonon frequency is depicted as an horizontal dotted line. Bottom panels: Plasmon (grey dashed lines), photon (magenta solid lines), and phonon (green dotted lines) admixtures of the antisymmetric lower polariton W LP i,qσ=− (i = pl, pt, pn) as a function of q/q 0 for ω pl = 1.5ω 0 . The top and bottom sub-panels correspond to ν pn = 0 and ν pn = 0.5ω 0 , respectively. regime where the symmetric SPP is mostly photon-like with w 0 qσ=+ ∼ qc. Interestingly, the resonance between the antisymmetric SPP and the phonon mode is now shifted to a large wave vector q r ≈ 220q 0 , where the antisymmetric SPP exhibits a pure plasmonic character.
For ν pn = 0, the antisymmetric LP (thick red solid line) and UP (thick blue solid line) are splitted in the vicinity of the resonance q r ≈ 220q 0 , while the symmetric polaritons (thick colored dashed lines) do not feature any anticrossing behavior for q ∼ q 0 . Similarly as in a), the LPs σ = ± can be seen as dressed phonon modes with frequencies red-shifted from ω 0 for large q due to their plasmon admixtures. Here, however, these dressed phonons can propagate with very large wave vectors comparable to the electronic Fermi wave vector.

Electron-phonon scattering
We now show that the ability to tune the energy and wave vector of the dressed phonons can affect the electron-phonon scattering in the crystal. The coupling Hamiltonian between individual electrons (electronic dark modes) and phonons reads The fermionic operators c K and c † K annihilate and create an electron with wave vector K and energy ξ K (K ≡ |K|) relative to the Fermi energy. For simplicity, we assume a 3D spherical , with m the electron effective mass and K F the Fermi wave vector, and that the coupling constant M does not depend on Q as shown theoretically for intramolecular phonons in certain crystals. 57 Electron-phonon interactions are usually characterised by the dimensionless coupling parameter λ, which quantifies the electron mass renormalization due to the coupling to phonons. 50 At zero temperature, λ is defined as where stands for real part, is the 3D electron density of states at the Fermi level, and ∂ ω Σ K (ω)| ω=0 denotes the frequency derivative of the retarded electron self-energy Σ K (ω) evaluated at ω = 0. An equation of motion analysis 58 of the electron where In the presence of phonon-photon coupling (ν pn = 0), the phonon dynamics is governed by the Hamiltonian H pol , which includes the coupling of phonons to plasmons and photons. It is therefore convenient to express the 3D phonon operators B Qα and B † Qα in terms of the 3D hybrid modes Π Qαj and Π † Qαj , and then project the latter onto the quasi-2D bright and dark modes such that the electron-phonon Hamiltonian Eq. (4) takes the form H el−pn = el−pn . The contribution of the bright modes reads el−pn of the dark modes that do not interact with photons is given in the supplemental material. Here, η j = χ j is associated to the hybrid modes weights of the phonons. The electron self-energy due to the interaction with bright modes reads where P qσj (ω) = −i dτ e iωτ Υ qσj (τ )Υ −qσj (0) and Υ qσj = π qσj + π † −qσj . As detailed in the supplemental material, we now use the equation of motion theory to calculate the GF since ω 1 = ω 0 (and ω 2 = ω pl ). In this case, η 1 = 1, η 2 = 0, and only the pure phonon term j = 1 therefore contributes to Σ (B) K . We then calculate P qσj (ω) and the self-energy Eq. (6) for ν pn = 0, and decompose the electron-phonon coupling parameter into its bright and dark mode contributions: λ = λ (B) + λ (D) for ν pn = 0, and λ 0 = λ can not solely lead to a modification of the electron-phonon scattering. Finally, the relative enhancement of the electron-phonon coupling parameter reads (see supplemental material) where the function ϕ qσ is given in the supplemental material. The latter describes the renormalization of the phonon energy due to the coupling to plasmons and photons, and depends on the polariton frequencies w qσζ . It is represented on Fig. 3 a) together with the frequency dispersion of the LPs w qσLP /ω 0 corresponding to the red solid and dashed lines on Fig. 2 b). First, we find that the contributions of the LPs to ϕ qσ largely dominate over the contributions of the other polaritons. Secondly, we observe that ϕ qσ is anticorrelated with the dispersion of the LPs, whenever the latter exhibit a finite phonon weight. For instance, ϕ qσ=− reaches large values in the region q q 0 where the frequency of the antisymmetric LP (∼ 35% phonon for ν pn = 0.5ω 0 ) is far away from the bare phonon frequency. In contrast, while w qσ=+LP ω 0 in this region, ϕ qσ=+ ≈ 1 since the symmetric LP is fully photon-like for q q 0 . Note that for ν pn /ω 0 → 0, the LPs for q > q r are fully phonon-like with w qσLP → ω 0 , and one finds that ϕ qσ → 1 does not contribute to ∆λ λ 0 in this region [see Eq. (7)]. Similarly, ϕ qσ does not contribute to ∆λ λ 0 for q < q r because of the vanishing phonon weight of the LPs for ν pn /ω 0 → 0.  Figure 3: a) Function ϕ qσ versus q/q 0 for ν pn /ω 0 = 0.5, q 0 = 0.01, and ω pl = 1.5ω 0 . Inset: Normalized frequency dispersion of the LPs w qσLP /ω 0 corresponding to the red solid and dashed lines on Fig. 2 b). The contributions of symmetric (σ = +) and antisymmetric (σ = −) modes are depicted as dashed and solid lines, respectively. b) Relative enhancement of the electron-phonon coupling parameter ∆λ λ 0 given by Eq. (7) versus q 0 and ν pn /ω 0 , for ω pl = 1.5ω 0 and K F = 10 3 q 0 .
The relative enhancement of the electron-phonon coupling ∆λ λ 0 is represented on Fig. 3 b) as a function of q 0 and ν pn /ω 0 . We find that it can reach large values 10% for ν pn 0.4ω 0 and sufficiently small (q 0 ∼ 0.01). The fact that ∆λ λ 0 decreases when increasing the crystal thickness can be understood by realizing that the function |f ασ (Q)| 2 entering Eq. (7) scales as ∼ 1/(γ cr ), as shown in the supplemental material. This means that the dressing of 3D phonons by evanescent waves decaying exponentially at each interface is large for thin enough crystals, when surface effects become important. Nevertheless, when the crystal thickness becomes too small, namely when K F 1, the description of the matter excitations by 3D fields breaks down and quantum confinement effects have to be taken into account in the model. The parameters of Fig. 3 b) are = 10nm for ω 0 = 0.2eV, and K F = 1nm −1 . We find that the enhancement of the electron-phonon coupling strongly increases with the ratio q 0 /K F , i.e. when the mismatch between the typical photonic and electronic wave vectors is reduced. This can occur in materials with larger phonon frequency and/or lower electron density. In the latter case, however, decreasing K F requires to consider larger for the 3D description of the matter fields to be valid, which in turn leads to a reduction of ∆λ λ 0 . For instance, bilayer graphene with ∼ 0.7nm 59 exhibits an infrared-active phonon mode at ω 0 = 0.2eV, 40 and has been recently shown to feature superconductivity at low carrier densities ≈ 2 × 10 12 cm −2 when the two sheets of graphene are twisted relative to each other by a small angle. 48 In this regime, we find that ω pl ≈ ω 0 , q 0 ≈ 7×10 −4 , K F ≈ 1.6×10 3 q 0 , and the large phonon-photon coupling strengths ν pn ≈ 0.25ω 0 40 allow to reach the ultrastrong coupling regime. Using these parameters, we find that the contribution of surface plasmonphonon polaritons to the electron-phonon scattering is expected to reach a very large value ∆λ λ 0 ∼ 2.7. However, the small thickness of bilayer graphene such that K F ∼ 1 in this weaklydoped regime implies that quantum confinement effects become important. Furthemore, the failure of the 3D model assuming ω pl → const as q → 0 can be seen in the dispersion of the optical plasmon mode in bilayer graphene, 53 which is known to be ω pl ∼ √ q as q → 0.

Conclusion
In summary, we have carried out a non-perturbative quantum theory of intrinsic, lossless

Supplemental material
In this supplemental material, we provide details of the calculations sketched in the main text.
The first section is dedicated to the derivation and diagonalization of the matter Hamiltonian leading to the hybrid plasmon-phonon modes. In the second section, we derive the photon and light-matter coupling Hamiltonians, and explain the self-consistent method used to find the surface polariton modes of the system. Further details concerning the theoretical foundations of the model can be found in Ref. [55]. The plasmon, phonon, and photon admixtures of the polaritons are provided at the end of the section. In the third section, we derive the coupling Hamiltonian between electronic dark modes and intramolecular optical phonons, and show in detail how the surface polariton modes can lead to an enhancement of the electron-phonon scattering in the crystal. In the last section, we study an elementary plasmonic structure composed of a single metal-dielectric interface, where phonons in the dielectric region strongly interact with SPPs. We use an effective model similar to that of Ref. [24], which is based on the SPP quantization scheme detailed in Refs. [60][61][62], and show that unphysical behaviors appear in the ultrastrong coupling (USC) regime.

Diagonalization of the matter Hamiltonian
We first explain the diagonalization procedure of H mat in detail. The latter involves both plasmonic and phononic degrees of freedom through the total polarization field P = P pl +P pn , where denotes the plasmon polarization field and P pn the phonon polarization. In the case of a 3D free electron gas, the electron-hole excitation operator P † KQ = c † K+Q c K is written in terms of the fermionic operators c K and c † K , which annihilates and creates an electron with 3D wave vector K, respectively. For long-wavelength excitations with Q K F , the Hilbert space is restricted to the RPA subspace spanned by the states {|F , P † KQ |F }, where |F = K<K F c † K |0 is the electron ground state (Fermi sea of the free electron gas), K = |K| the wave vector modulus, |0 the vacuum state, and K F the Fermi wave vector. In this case, it can be shown 55 that P † KQ and its hermitian conjugate are bosonic operators satisfying the commutation relation [P KQ , P † K Q ] = δ K,K δ Q,Q . The electronic dipole moment reads where n K denotes the Fermi occupation number at zero temperature, namely n K = 1 for K < K F and n K = 0 for K > K F . The electron effective charge and mass are denoted as e and m, respectively. Denoting the lattice site positions as R i , where i ∈ [1, n] with n the number of unit cell in the crystal, the phonon polarization field can be written as Here with a 3 the volume of a unit cell. The matter Hamiltonian can be decomposed as H mat = H pl +H pn +H P 2 . The contribution H pl = K,Q ∆ξ K,Q P † KQ P KQ is the effective Hamiltonian providing the energy ∆ξ K,Q = ξ K+Q − ξ K of the electron-hole pairs created across the Fermi sea. The electron energy (relative to the Fermi energy) is denoted as ξ K . For simplicity, we assume a 3D spherical Fermi surface providing ξ K = (K 2 − K 2 F )/(2m). The free phonon Hamiltonian reads H pn = Q,α ω 0 B † Qα B Qα , and the interaction term H P 2 is proportional to the square polarization field: This term contains the contributions of both plasmons and phonons ∝ P 2 pl and ∝ P 2 pn , respectively, as well as a direct plasmon-phonon interaction term ∝ P pl · P pn . Replacing Eqs. (8) and (10) into Eq. (11) and performing the spatial integration, the matter Hamiltonian is derived as where we have introduced M 0 a 3 the ion plasma frequency. The Hamiltonian Eq. (12) can be diagonalized using the Hopfield-Bogoliubov method, 63 by introducing collective eigenmodes in the form The condition for these modes to diagonalize the Hamiltonian Eq. (12) is then [Π Qj , H mat ] = ω j Π Qj . Computing this commutator and introducing the vector we obtain the transformation coefficients entering Eq. (13): represents the transverse dielectric function of the crystal. To obtain this result, we have used the definition of the electronic plasma frequency ω pl = ρe 2 m 0 with ρ the electron density, α u α · u Q,p = 1, as well as the identity in the so-called dynamical long-wavelength limit Q Fermi velocity). Note that in this regime, one can use the relation to put Eq. (17) in the form of the well-known Lindhard function, which describes the longitudinal response of a free electron gas. 64 The two eigenvalues ω j can be determined by solving the equation cr ( ω j ) = 0, which provides the solutions given by Eq. (1) of the main text. As detailed in Ref. [55], the vector Z Qj can be decomposed as Z Qj = N j T j u Qp when solving for the transverse modes, where T j are determined by computing the residues of the inverse dielectric function 1 cr (ω) and are derived as The constant N j can be found by using the normalization condition stemming from the commutation relation [Π Qj , Π † Q j ] = δ Q,Q δ j,j . Using Eq. (15) with Z Qj = N j T j u Qp , as well as the identity which can be derived similarly as Eq. (17) in the dynamical long-wavelength limit, we obtain Finally, we use Eqs. (8), (10), (13), (15), the decomposition Z Qj = N j T j u Qp , as well as the to write the total (transverse) polarization field as Since the crystal is assumed to be isotropic, one can choose the transverse polarization vectors u Qp to coincide with u α (α = z,//), and using Eqs. (18) and (19) Denoting by u ⊥ the in-plane unit vector perpendicular to both u z and u // , and θ the heaviside function, the mode functions read Replacing these expressions into H pt = 1 2 0 dR D 2 (R)+ 1 2 0 c 2 dR H 2 (R) and performing the integrations, the free photon Hamiltonian is derived as with the overlap matrix elements after spatial integration. The function F ± (Q) stems from the overlap between the displacement and the polarization fields and reads F ± (Q) = (F 2 (Q) ± F 1 (Q)) / √ 2 with We now introduce the quasi-2D "bright" modes π † qσj = qz,α f ασ (Q)Π † Qjα and its hermitian conjugate, which consist of linear superpositions of the 3D plasmon-phonon hybrid modes.
The function f ασ (Q) reads The normalization constant N qσ is determined by imposing the bosonic commutation rela-tions [π qσj , π † q σj ] = δ q,q , which provides 1/N 2 q± = qz q 2 |F ± (q z )| 2 + γ 2 cr |F ∓ (q z )| 2 . Transforming the summation into an integral qz → /(2π) +∞ −∞ dq z , and using the results and we finally obtain Note that the existence of a symmetry plane at z = /2 implies that the symmetric and antisymmetric bright modes commute, i.e. [π qσj , π † q σ j ] = δ q,q δ σ,σ δ j,j , which can be checked using Eqs. (23), (24), and (25). On the other hand, the commutation relation with respect to j follows directly from the diagonalization of the matter Hamiltonian presented before. One can then use Eq. (26) combined with the definition of the bright modes to write the light-matter coupling Hamiltonian Eq. (22) in the form given in the main text. Using the unitary transformation π † qσj = qz,α f ασ (Q)Π † Qjα and Γ † qσjs = qz,α g ασs (Q)Π † Qjα , with [Γ qσjs , Γ † q σ j s ] = δ q,q δ σ,σ δ j,j δ s,s , [π qσj , Γ † q σ j s ] = [π † qσj , Γ † q σ j s ] = 0, and the identities qz,α g ασs (Q)g * ασ s (Q) = δ s,s δ σ,σ qz,α the matter Hamiltonian can be decomposed as H mat = q,σ,j ω j π † qσj π qσj + H dark , where s=1 ω j Γ † qσjs Γ qσjs denotes the contribution of the quasi-2D dark modes which do not interact with light. Leaving this contribution aside, the polariton Hamiltonian H pol of the main text can be diagonalized numerically using a Hopfield-Bogoliubov transformation. 63 We introduce the polariton eigenmodes of the system in the form which satisfy the eigenvalue equation [P qσζ , H] = w qσζ P qσζ . Computing this commutator, one can show that the polariton frequencies w qσζ in each subspace (q, σ) are given by the 3 positive eigenvalues of the matrix Out of these three eigenvalues, only the lowest two correspond to surface modes as located below the light cone. These two surface modes are denoted as lower and upper polaritons with frequencies w qσLP and w qσUP , respectively. In order to compute these frequencies, we use a self-consistent algorithm which consists in starting with a given frequency w qσLP , then determine the penetration depths in each media from the Helmholtz equation n (w qσLP )w 2 qσLP /c 2 = q 2 −γ 2 n with d = 1 and cr given by Eq. (16), and use the values of γ n to compute the parameters entering the matrix Eq. (28). The latter is diagonalized numerically allowing to determine the new w qσLP , and the whole process is repeated until convergence is reached. This method is applied independently for the symmetric and the antisymmetric modes σ = ±.
The photon, phonon, and plasmon weights of the polaritons modes are directly related to the transformation coefficients entering Eq. (27). From the polariton normalization the photon weight is defined as W ζ pt,qσ = −2iC qσ X qσζ Y * qσζ , and we now need to compute the phonon and plasmon weights of the hybrid plasmon-phonon modes j = 1, 2. This can be done by realizing that the frequencies ω 1 and ω 2 given by Eq. (1) are the eigenvalues of the effective plasmon-phonon coupling Hamiltonian H eff = ω pl P † P + ω 0 B † B + Ω P + P † B + B † , where P and B denote respectively the effective plasmon and phonon bosonic operators, Ω = (ν pn /2) ω pl / ω 0 is the effective coupling strength, and ω 0 = ω 2 0 + ν 2 pn the longitudinal phonon frequency. We introduce the normal modes P 0 satisfy the eigenvalue equation [P 0 j , H eff ] = ω j P 0 j . The coefficients U j , U j , V j , and V j can then be computed similarly as before by diagonalizing the associated Hopfield matrix. Using the normalization condition |U j | 2 − | U j | 2 + |V j | 2 − | V j | 2 = 1, the phonon and plasmon weights of the polariton ζ are defined as and satisfy the sum rule i W ζ i,qσ = 1 with i = pt, pl, pn.

Electron-phonon scattering
In this section, we derive the expression of the coupling Hamiltonian H el−pn between electronic dark modes and intramolecular phonons, 57,65,66 and show how the electron-phonon coupling parameter λ is affected by the phonon dressing. The potential generated by the vibrating ions in the crystal is expanded in the vicinity of the equilibrium is proportional to the phonon polarization field given by Eq. (10), and reads

Introducing the fermion field
a Wannier function localized on site i, the electron-phonon coupling Hamiltonian reads where the bare phonon contribution is contained in the polariton Hamiltonian derived in the previous section. Using these expressions and retaining only the diagonal terms evaluated at the same position (intramolecular vibrations), H el−pn takes the form with the coupling constant M = eV 00 N V 2 M ω 0 a 3 independent on the wave vector, and V 00 = dR |φ(R)| 2 ∇ V (R) · u α . Electron-phonon interactions are often characterised by the dimensionless coupling parameter λ, which quantifies the electron mass renormalization due to the coupling to phonons. 50,67 At zero temperature, λ reads where stands for real part, is the 3D electron density of states at the Fermi level, and ∂ ω Σ K (ω)| ω=0 denotes the frequency derivative of the retarded electron self-energy Σ K (ω) evaluated at ω = 0. An equation of motion analysis 58 of the electron Green's function (GF) G K (τ ) = −i c K (τ )c † K (0) allows to write the electron self-energy as where B Qα (ω) = −i dτ e iωτ B Qα (τ )B −Qα (0) is the phonon GF written in the frequency domain, and B Qα = B Qα + B † −Qα . In the absence of phonon-photon coupling (ν pn = 0), there is no hybridization between phonons and plasmons, nor is there coupling to photons. Phonons therefore enter the Hamiltonian H pol only in the free contribution H pn = Q,α ω 0 B † Qα B Qα . In this case, the equation of motion analysis simply provides B Qα (ω) = 2ω 0 /(ω 2 − ω 2 0 ). Using this expression together with the non-interacting electron GF G 0 Eq. (31), the retarded electron self-energy is derived as 50 where 0 + denotes a vanishing positive number. As detailed in the chapter 6 of Ref. [50], we keep only the terms ∝ n K−Q in Eq. (32), and calculate the quantity ∂ ω Σ K (ω)| ω=0 . We then use an integration by parts and consider only the term ∝ dQ We finally obtain the electron-phonon coupling parameter Eqs. (30) for ν pn = 0: In the presence of phonon-photon coupling (ν pn = 0), the phonon dynamics is governed by the Hamiltonian H pol , which includes the coupling of phonons to plasmons and photons. Using Eqs. (13), (15), (18), and (19), we express the 3D phonon operators B Qα and B † Qα in terms of the 3D hybrid modes Π Qαj and Π † Qαj , and then project the latter onto the quasi-2D bright and dark modes defined in the previous section such that the electron-phonon Hamiltonian el−pn . The contribution of the bright modes reads ) is associated to the hybrid modes weights of the phonons. Similarly, the contribution of the dark modes is Using the Hamiltonians H el−pn , the electron self-energies due to the interaction with bright and dark modes are derived as Here, P qσj (ω) and G qσjs (ω) denote the bright and dark mode GFs, respectively, defined as P qσj (ω) = −i dτ e iωτ Υ qσj (τ )Υ −qσj (0) and G qσjs (ω) = −i dτ e iωτ Ξ qσjs (τ )Ξ −qσjs (0) , with Υ qσj = π qσj +π † −qσj and Ξ qσjs = Γ qσjs +Γ † −qσjs . We now decompose the electron-phonon coupling parameter into its bright and dark mode contributions, which are respectively associated to the self-energies Eqs. (33) and (34): λ = λ (B) + λ (D) for ν pn = 0, and λ 0 = in the case ν pn = 0. For ν pn = 0, the bright and dark mode GFs read P qσj (ω) = G qσjs (ω) = 2 ω j /(ω 2 − ω 2 j ), with ω 1 = ω 0 , ω 2 = ω pl , η 1 = 1, and η 2 = 0. Using these GFs together with G 0 K (ω) = 1/(ω − ξ K ) in Eqs. (33) and (34), we proceed as before and find For ν pn = 0, the dark modes become hybrid plasmon-phonon modes, but do not interact with photons. Still using the same procedure with the GF G qσjs (ω) = 2 ω j /(ω 2 − ω 2 j ) in Eq. (34) and G 0 K (ω) = 1/(ω − ξ K ), we obtain Interestingly, the sum rule We now use the equation of motion theory 58,68 to calculate the bright mode GF P qσj (ω) for ν pn = 0: We start from the Heisenberg equation of motion and calculate the double commutator. A vanishing source term q,σ,j J qσj Υ −qσj is then added to H pol , and we calculate the functional derivative of the obtained equation of motion with respect to J qσj (0). The bright mode GF can be written as where · · · 0 denotes the expectation value in the non-interacting ground state of H pol (without the light-matter coupling term H mat−pt ). Introducing the photon GF D qσ (τ ) = δ Dqσ(τ ) δJ qσj (0) , we obtain the following set of equations in the frequency domain: where D 0 qσ (ω) = 2Wqσ ω 2 −W 2 qσ and P 0 qσj (ω) = 2 ω j /(ω 2 − ω 2 j ) denote the non-interacting photon and bright mode GFs, and W qσ = 2cC qσ α qσ β qσ plays the role of the photon frequency.
Solving for P qσj (ω), we obtain with the photon self-energy S jj qσ (ω) = 2cC 2 qσ βqσΩ qσj Ω qσj . Note that the polariton frequencies w qσζ correspond to the poles of P qσj (ω), and can be determined by solving the associated sixth-order polynomial equation. Proceeding as before, we use the GF Eq. (38) together with G 0 K (ω) = 1/(ω − ξ K ) in Eq. (33), as well as Eqs. (36) and (37), and obtain the relative enhancement of the electron-phonon coupling parameter: where the function ϕ qσζ describes the renormalization of the phonon energy due to the coupling to plasmons and photons, and depends on the polariton frequencies as One can then transform the summation over Q into an integral, and write the relative enhancement in the simple form ∆λ

Effective quantum model in the USC regime
In this section, we show that an effective quantum model similar to that of Ref. [24] exhibits unphysical behaviors in the ultrastrong coupling regime. We consider a single metal-dielectric interface featuring a plasmon mode with plasma frequency ω m in the metal region, and a collection of vibrating ions with frequency ω 0 forming a lattice embedded in the dielectric region of permittivity 1 (e.g. air). Phonons are polarized in the directions u z and u // , and both the metal and dielectric regions extend to infinity on each side. In the effective quantum model of Ref. [24], the Power-Zienau-Woolley Hamiltonian of the system reads H pol = H spp + H pn−spp + H pn + H P 2 , where H spp is the effective SPP Hamiltonian, H pn−spp is the SPP-phonon coupling, and H pn = Q,α ω 0 B † Qα B Qα is the phonon Hamiltonian. The term proportional to the square polarization reads H P 2 = 1 2 0 dR P 2 pn (R), where the phonon polarization field P pn is given by Eq. (10). One can show that the overall phonon contribution to the Hamiltonian can be put in the diagonal form H pn + H P 2 = ω 0 Q,α B † Qα B Qα , with the transformation The longitudinal phonon frequency is defined as ω 0 = ω 2 0 + ν 2 pn with ν pn the ion plasma frequency. The effective SPP Hamiltonian reads In contrast to the model presented in the previous sections, here, the displacement and magnetic fields D and H contain the contribution of the plasmon polarization and do not involve pure photonic degrees of freedom. The displacement field can be written as with S an arbitrary (large) surface in the plane, and a q and a † q the annihilation and creation operators of a SPP with frequency w 0 q , which satisfy the bosonic commutation relations [a q , a † q ] = δ q,q . The mode functions read , with m (ω) = 1 − (ω m /ω) 2 the Drude permittivity of the (lossless) metal. According to the SPP quantization scheme detailed in Ref. [62], the effective quantization length L eff (ω) is defined as With these definitions, one can show 60-62 that the effective SPP Hamiltonian takes the usual form H spp = q w 0 q a † q a q . Expressing the phonon polarization field Eq. (10) with the vacuum Rabi frequency and can be diagonalized using a Hopfield-Bogoliubov transformation. The two resulting polaritons modes are represented on Fig. 4 for ν pn = 0.01ω 0 (dotted lines), ν pn = 0.5ω 0 (dashed lines), and ν pn = ω 0 (solid lines). The lower ζ = LP and upper ζ = UP polaritons with frequency w qζ are depicted as red and blue lines, respectively. While the limiting behaviors of the LP w qLP ∼ qc for q → 0 and the UP w qUP → ω m / √ 2 for q → ∞ are properly described by the model, w qUP for q → 0 and w qLP for q → ∞ feature unphysical behaviors in the ultrastrong coupling regime, i.e. when ν pn is a non-negligible fraction of ω 0 . Indeed, w qLP exhibits irregularities in the regime of large detunings q q 0 (see inset), and becomes larger than ω 0 while increasing the coupling strength ν pn . We conclude that the effective quantum model based on quantized SPPs is well suited as long as the coupling strength does not exceed a few percent of the associated transition frequency, and therefore fails in the ultrastrong coupling regime. : Log-log scale: Normalized polariton frequencies w qζ /ω 0 versus in-plane wave vector q/q 0 (q 0 = ω 0 /c), for ν pn = 0.01ω 0 (dotted lines), ν pn = 0.5ω 0 (dashed lines), and ν pn = ω 0 (solid lines). The lower ζ = LP and upper ζ = UP polaritons are depicted as red and blue lines. The region delimited by the rectangle is magnified in the inset. The plasma frequency is chosen as ω m = 45ω 0 , which corresponds to ω m = 9eV (gold or silver) for ω 0 = 200meV (mid-infrared optical phonons).