Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity
Michael A. Sentef Michael Ruggenthaler Angel Rubio






[]michael.sentef@mpsd.mpg.de
Max Planck Institute for the Structure and Dynamics of Matter,
Center for Free Electron Laser Science, 22761 Hamburg, Germany


Max Planck Institute for the Structure and Dynamics of Matter,
Center for Free Electron Laser Science, 22761 Hamburg, Germany


Max Planck Institute for the Structure and Dynamics of Matter,
Center for Free Electron Laser Science, 22761 Hamburg, Germany
 Nano-Bio Spectroscopy Group, Universidad del País Vasco, , 20018 San Sebastián, Spain 
Center for Computational Quantum Physics (CCQ), The Flatiron Institute, 162 Fifth Avenue, New York NY 10010








Laser control of solids was so far mainly discussed in the context of strong classical nonlinear light-matter coupling in a pump-probe framework. Here we propose a quantum-electrodynamical setting to address the coupling of a low-dimensional quantum material to quantized electromagnetic fields in quantum cavities. Using a protoypical model system describing FeSe/SrTiO_3 with electron-phonon long-range forward scattering, we study how the formation of phonon polaritons at the 2D interface of the material modifies effective couplings and superconducting properties in a Migdal-Eliashberg simulation. We find that through highly polarizable dipolar phonons, large cavity-enhanced electron-phonon couplings are possible but superconductivity is not enhanced for the forward-scattering pairing mechanism due to the interplay between coupling enhancement and mode softening. An analysis of critical temperature dependencies on couplings and mode frequencies suggests that that cavity-enhanced superconductivity is possible for more conventional short-range pairing mechanisms. Our results demonstrate that quantum cavities enable the engineering of fundamental couplings in solids paving the way to unprecedented control of material properties.



Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity


    A. Rubio
    December 30, 2023
===========================================================================================================================








§ INTRODUCTION

Strong coupling and manipulation of matter with photons in quantum-electrodynamical (QED) environments is becoming a major research focus across many disciplines. Among the topics with large potential are the creation of exciton-polariton condensates <cit.>, polaritonic chemistry <cit.> and transport <cit.>, quantum nanoplasmonics <cit.>, light-induced topology <cit.> and magnetism in 2D materials <cit.>, and novel spectroscopies <cit.>. In condensed matter, the search for control knobs that allow to design properties of quantum materials is an ongoing broad research effort <cit.>. One possible route is to employ the nonequilibrium dynamics and coherent manipulation of quantum many-body systems with ultrashort laser pulses <cit.>. However, in these cases “classical” light was typically used.  Here we propose a new route towards manipulating microscopic couplings in solids and inducing ordered phases especially at interfaces and in two-dimensional materials. 

The discovery of enhanced superconductivity in monolayer FeSe on SrTiO_3 <cit.> and its possible relation to a cross-interfacial electron-phonon coupling <cit.> has stimulated considerable interest with an ongoing open debate <cit.>. Irrespective of the outcome of this debate, the interfacial phonon mode under consideration is of particular interest for light-control purposes as it has a dipole moment implying bilinear phonon-photon coupling, while at the same time the phonon also couples bilinearly to in-plane FeSe electrons, with a vertex that is strongly peaked for small momentum transfers known as forward scattering. This combination of features is due to the high degree of anisotropy owing to the interfacial structure. Here we employ a prototypical model system, related to FeSe/SrTiO_3, for such extreme forward scattering to investigate how photon-phonon coupling in cavities can affect electron-phonon coupling and phonon-mediated superconductivity.




§ RESULTS




 §.§ Setup: Two-dimensional material inside a cavity


In Figure <ref>A we show the setup for a two-dimensional material inside a QED cavity environment with perfectly reflecting mirrors. The mirrors confine the photon modes inside the cavity and can lead to strong light-matter coupling even when only vacuum of the electromagnetic field is considered <cit.>. Specifically, we propose a layered structure of a 2D material (e.g., monolayer FeSe) on a dielectric substrate with large dielectric constant (e.g., SrTiO_3) that further helps confine the cavity photon modes of interest. 

For the particular example of FeSe/SrTiO_3, the effect of the cavity is to couple the electromagnetic field of the photons polarized along the z direction, perpendicular to the interfacial plane, to a cross-interfacial phonon mode. Importantly, here we go beyond the often-employed rotating-wave and dipole approximations for the light-matter interaction and use full minimal dipolar coupling including the J · A and A^2 terms (see Section B of the Supplementary Materials), which makes the theory manifestly gauge-invariant and avoids unphysical divergences. The phonon has a dipole moment along z that involves motion of the O and Ti ions in the topmost layer of SrTiO_3, spatially very close to the FeSe monolayer. Specifically, one quasi-dispersionless optical Fuchs-Kliewer phonon at 92 meV <cit.> was identified as the most relevant phonon mode that strongly couples to the FeSe electrons both in angle-resolved photoemission <cit.> and high-resolution electron energy loss spectroscopies <cit.>. The influence of screening on this mode is not settled yet, in particular when it comes to phonon linewidths <cit.>. However, the experimental evidence for its influence on electronic properties <cit.> is definitely present suggesting to use this mode to build a simplified model Hamiltonian to address the impact of reaching strong light-matter coupling on the superconducting behavior of the material. We specifically use a single-band model for the electrons in two spatial dimensions in a partially filled band with filling n=0.07 per spin, as previously used to model the relevant electronic structure fitting angle-resolved photoemission data <cit.>. A bilinear electron-phonon scattering is introduced by a coupling vertex g()=g_0 exp(-||/q_0) that is strongly peaked near momentum =0 with a coupling range q_0. The coupling strength g_0 is adjusted to keep a total dimensionless coupling strength λ≈ 0.18 independent of q_0, where λ is determined from the effective electronic mass renormalization m^*/m = 1+λ in the metallic normal state above the superconducting critical temperature in absence of the cavity coupling. This conservative choice of λ is for instance below the value of 0.25 that was given in Ref. zhang_role_2016.

Through phonon-photon coupling we study phonon-polariton formation in this setting. In Figure <ref>B we show schematically the resulting polariton branches that stem from a gauge-invariant coupling involving both J · A and A^2 terms, where J is the current of phononic dipoles associated to an infrared-active phonon mode, and A the electromagnetic gauge field of the photons. The relevant effective coupling strength between photons and phonons is given by the phononic plasma frequency ω_P = √(4 π e^2/M ν_0, 2D L_z), with M the reduced mass of the phonon (see Section B of the Supplementary Materials). For the 2D system in the cavity the plasma frequency is controlled by the length of the vacuum inside the cavity in z direction, L_z, and the 2D unit cell area ν_0, 2D=L_x L_y / N_x N_y, with L_i, N_i the length and number of unit cells of the system in i direction, respectively. The plasma frequency sets the splitting between the upper and lower polariton branches, reminiscent of the LO-TO splitting in bulk semiconductors. Obviously this splitting is only relevant at very small momenta q since the photon energies become large compared to the phonon frequency quickly as q increases due to the large magnitude of the speed of light. 

The formation of phonon polaritons leads to a redistribution of the electron-phonon coupling vertex into the two polariton branches. In the following, we refer to this coupling between electrons and phonon polaritons as “electron-phonon coupling”, since the coupling originates from electron-phonon coupling in the free-space setting without cavity, and direct electron-photon coupling is not relevant in our setup. In Figure <ref>C we plot the squares of the coupling vertices between electrons and the respective polaritons as a function of q/k_F, where k_F is the Fermi momentum. A realistic value of the coupling range for FeSe/SrTiO_3 was estimated as q_0/k_F ≈ 0.1, as needed to create replica bands in angle-resolved photoemission that duplicate primary band features without significant momentum smearing <cit.>. In a microscopic model, this value depends on the distance h_0 between the topmost TiO_2 layer and the FeSe monolayer as well as the anisotropy of in-plane and perpendicular dielectric constants via q_0^-1 = h_0 √(ϵ_∥/ϵ_), with realistic estimates ϵ_∥/ϵ_≈ 100 and 1/(h_0 k_F) ≈ 1. This coupling range is larger than the momentum at which photon and phonon branches cross and mix most strongly in the polariton formation process. This implies that the modification of electron-phonon coupling due to the cavity only happens at very small momenta typically smaller than q_0/k_F. Thus to investigate how the degree of forward scattering influences the way in which cavity coupling is able to modify the electronic properties, we employ different values for q_0/k_F below, envisioning that cavity effects are enhanced when q_0/k_F becomes smaller, which would in practice be achieved by making the dielectric-constant anisotropy ratio larger. In Table <ref> we summarize the relevant parameter values of the bare material used in our simulations.





 §.§ Cavity-enhanced electron-phonon interaction

The critical question to answer here is how the redistribution of the coupling vertex to the upper and lower polariton branches affect the electronic properties. We investigate this by a diagrammatic approach employing Matsubara Green's functions. We adopt the same approximations used in Ref. rademaker_enhanced_2016 and compute the self-consistent Migdal-Eliashberg diagram with dressed electronic Green's function in Nambu space, allowing us to take into account superconducting order. The central quantity is the electronic self-energy Σ̂(,iω_n) =
    iω_n[1-Z(,iω_n)]τ̂_0 + χ(,iω_n)τ̂_3 +
    ϕ(,iω_n)τ̂_1,
written in terms of the Pauli matrices τ̂_i, the effective mass renormalization Z(,iω_n), the band dispersion renormalization χ(,iω_n), and the anomalous
self-energy ϕ(,iω_n), which vanishes in the normal state. 

We first investigate the effect of the cavity on the effective electron-phonon coupling λ itself. This is of interest independently of superconductivity to be discussed below, as the electron-phonon coupling affects also many other properties of materials, such as the conductivity, structural phase transitions, or superconductivity in standard BCS superconductors. In particular it plays a pivotal role for THz-driven nonequilibrium phases of materials. In Figure <ref> we show how cavity coupling modifies the temperature-dependent quasiparticle mass renormalization obtained from the normal self-energy for the different coupling ranges, realistic q_0/k_F=0.105 (Figure <ref>A), reduced q_0/k_F=0.053 (Figure <ref>B), and very small q_0/k_F=0.021 (Figure <ref>C). The first observation is that independent of the cavity λ shows a strong temperature dependence with a peak around T_C, decreasing both towards higher temperatures and towards lower temperatures deep inside the ordered phase. The former is readily understood as a usual temperature effect when at high temperature the system becomes more and more classical and less correlated. The latter is understood by considering the fact that correlation effects are reduced deep in the ordered phase when quantum fluctuations lose their importance and a quasi-classical mean-field description can be adopted. Importantly, λ is enhanced by the cavity at all temperatures. The cavity effects are more pronounced as ω_P increases for fixed q_0/k_F, and as q_0/k_F increases for fixed ω_P.



 §.§ Light-modified superconductivity

We now turn to the effect of the cavity on superconductivity. Naively one might expect that an enhanced λ leads to enhanced superconducting critical temperature T_C. However, the relation is nontrivial as also the effective polariton frequency is relevant for T_C. We will see in the following that, unfortunately, for our system the enhancement of λ is cancelled by a reduction in the effective frequency.

Figure <ref>A shows the resulting temperature-dependent superconducting order gap Δ≡ϕ(_F, iπ/β)/Z(_F, iπ/β) evaluated at the smallest Matsubara frequency and at a Fermi momentum _F ≈ (0.666/a, 0.666/a) along the Brillouin zone diagonal for a coupling range q_0/k_F = 0.105 representative of FeSe/SrTiO_3. Starting from a critical temperature T_C ≈ 63 K in the absence of the cavity (ω_P=0.0), we find a slight reduction of superconductivity as the cavity is introduced and its extension L_z in the z direction perpendicular to the 2D material is reduced, resulting in a nonzero ω_P∝ 1/√(L_z). For perhaps unrealistically large values ω_P = 5.0 (eV), a reduction of T_C on the order of 1 Kelvin is found in our simulations, which would likely require cavity sizes of a few lattice constants and might in practice be too small to achieve at the moment.

In order to investigate the effect of the forward-scattering coupling range, we look at the change of the superconducting order in the case of q_0/k_F = 0.053 that is reduced by a factor of two from the realistic value described above, see Figure <ref>B. In this case the polaritonic redistribution of the coupling is expected to be more effective as there is a better match between the coupling range and the polariton mixing. This is indeed observed in the superconducting order enhancement. Where a value of ω_P = 5.0 was needed in Figure <ref>A to obtain a visible modification of T_C, here a smaller value ω_P = 2.5 is sufficient to enhance T_C by ≈ 1K. Even larger ω_P lead to enhancements of order 5 %. Finally if we decrease the range by another factor of two, q_0/k_F=0.021, the modification is relatively strong with changes of more than 10 %, shifting T_C by up to 10 K (Figure <ref>C).



 §.§ Analysis of the influence of the cavity on superconductivity

In order to gain physical intuition into why the enhancement of λ is insufficient to enhance superconductivity, we take a look at the approximate equation for T_C derived by Rademaker et al. <cit.> in the extreme forward-scattering and weak-coupling limit:

    T_C    ≈λΩ/2 + 3 λ.

From this expression it becomes clear that the enhancement of λ has to be sufficiently strong compared to the suppression of Ω that happens concomitantly in our case. This should be contrasted with the standard expression for a momentum-independent coupling vertex in Bardeen-Cooper-Schrieffer (BCS) theory,
T_C, BCS≈ 1.13 Ωexp(-1/λ).
The quasi-linearity in λ in Eq. (<ref>) leads to relatively high T_C for moderate values of λ, but in the cavity also has the negative effect that the enhancement of T_C scales only linearly rather than exponentially with λ.




§ DISCUSSION

Unfortunately, the enhancement of λ predicted here does not lead to an enhancement of the superconducting critical temperature T_C in our chosen setting. This effect is explained by the linear scaling of the critical temperature with λ for the case of extreme forward scattering in contrast to the exponential scaling for momentum-independent coupling. However, for more conventional pairing mechanisms not geared towards forward scattering, the observed enhancement of λ could naturally lead to enhanced T_C. Moreover, our theory and the analytical estimates of T_C are valid only in the Migdal-Eliashberg regime of weak coupling, unrenormalized polaritons, and adiabaticity. A polaritonic enhancement of λ could still lead to enhancement of T_C even for the forward-scattering case at intermediate couplings, when feedback effects on the polaritons become important, and when nonadiabatic effects come into play. Similarly, interplay between polaritonic pairing and other pairing mechanisms such as spin or orbital fluctuations are subjects for future study. It is possible that in such cases our original motivation of this work, namely to enhance superconductivity in a cavity, might work out. 

In summary we propose to employ QED cavity settings to control polaritonically mediated effects in low-dimensional materials. In reality the size of the achieved effects will depend on the quality factor of the cavity, the degree to which our idealized boundary conditions are realized in practice, and on the required large coupling strengths that can actually be reached in real devices. Importantly, however, our above results are ground-state modifications that are still qualitatively valid even in dissipative systems <cit.>. Moreover for organic molecules in cavities the ultrastrong-coupling regime was even achieved in bad cavities with small quality factors <cit.>. Here we predict changes of T_C in a few percent range for few-percent changes of the electron-phonon coupling λ. 
Known examples of LO-TO splitting in bulk semiconductors such as GaP suggest typical ratios of ω_P/Ω of order 10 % <cit.>, an order of magnitude smaller than the ones employed in this work. However, we caution that these are very different materials from the ones employed here, and oxide dielectrics close to the ferroelectric phase transition, such as SrTiO_3, were suggested to have giant LO-TO splittings exceeding 50% of the TO frequency <cit.> due to enhanced Born effective charges placing them much closer to the values explored here. It remains to be answered how large realistic LO-TO splittings can become at interfaces. It will definitely be important to explore strategies for enhancing the plasma frequency by synthesizing samples using different substrates with strongly coupled polar phonons, and exploring interface and heterostructure engineering to optimize the dielectric environment.  

We note that a related idea of exciton-mediated superconducting pairing <cit.> in 2D heterostructures was introduced <cit.> and recently discussed in the context of transition-metal dichalcogenides <cit.>. These proposals require exciton-polariton condensates to exist in the first place, which then affect pairing in doped nearby layers via coupling of quasifree electrons to condensed exciton polaritons. By contrast, our present proposal does not rely on bosonic condensation but rather focusses on directly modifying the electron-phonon coupling through polariton formation in a cavity. For the example of FeSe/SrTiO_3, our proposal could help shed light on the above-mentioned debate about the role of the forward-scattering phonon for superconductivity. If the coupling of the phonon to electrons is unimportant, the polaritonic effects will not play a role, which could serve as a test for the influence of the phonon on the electronic properties. Similarly, it was recently suggested to use classical lasers in a pump-probe setting to study the forward-scattering nature of the phonon <cit.>. Ongoing work focuses on a realistic ab initio computation of cavity-enhanced couplings via dipolar phonons using the framework of quantum-electrodynamical density functional theory <cit.>.

Note added in revision. Upon revision of the manuscript we became aware of two related works, that discuss related ideas of modifying superconducting properties by electron-photon interactions in cavities <cit.>. 




§ MATERIALS AND METHODS


We employ a cavity quantum-electrodynamical setting with plane-wave mode expansion inside a cavity, with fixed-node boundary conditions for confined cavity photon modes along the z direction, and periodic boundary conditions in the extended 2D plane (see Section A of the Supplementary Materials). Specifically, we use the Migdal-Eliashberg approximation to the electronic self-energy to a coupled electron-polariton model Hamiltonian involving electron-phonon forward scattering and dipolar phonon-photon coupling.

The electron-polariton Hamiltonian has the form 

    H     = ∑_,σϵ^†_ c^†_,σc^†_,σ
        + 1/√(N)∑_,,σ,λ=± c^†_+,σ
        c^†_,σ (g^*_λ() α^†_-,λ + g_λ() α^†_,λ) 
        + ∑_q⃗, λ=±ω_λ(q⃗) α^†_q⃗,λα^_q⃗,λ,

with c^†_,σ (c^†_,σ) the electron creation (annihilation) operators at
wavevector  and spin σ, ϵ_ = -2t[cos(k_x a)+cos(k_y a)] - μ the electronic band dispersion measured relative to the chemical potential
μ which is adjusted to fix a band filling of 0.07 per spin. Furthermore N is the number of k points in the 2D Brillouin zone, and g_λ() is the polariton-momentum -dependent electron-polariton coupling to branch λ=±,

    g_+()    = i sin(θ_) √(ω_+(q⃗)/Ω)  g_0 exp(-||/q_0), 
    
    g_-()    = i cos(θ_) √(ω_-(q⃗)/Ω)  g_0 exp(-||/q_0),

with bosonic polariton creation (annihilation) operators α^_q⃗,λ (α^†_q⃗,λ) for the polaritons with energies

    ω_±()    = (1/2( ω_phot(q⃗)^2 + ω_P^2 + Ω^2 ±√((ω_phot(q⃗)^2 + ω_P^2 + Ω^2)^2 - 4 ω_phot(q⃗)^2 Ω^2
     )) )^1/2.

The unitary transformation from phonons and photons to polaritons is parametrized by

    arctan(θ_)    = ω_phot(q⃗)^2 + ω_P^2 - Ω^2 + √((ω_phot(q⃗)^2 + ω_P^2 + Ω^2)^2 - 4 ω_phot(q⃗)^2 Ω^2 )/2 Ωω_P.

Here the underlying bare energies are given by the electronic hopping t = 0.075eV <cit.>, the phonon frequency Ω = 92 meV <cit.>, the bare photon dispersion is ω_phot(q⃗) = c || with speed of light c, and we use a variable effective phononic plasma frequency ω_P throughout the main text. Further details can be found in Sections B and C of the Supplementary Materials.

The Migdal-Eliashberg electronic self-energy on the Matsubara frequency axis is given by

    Σ̂(,iω_n) =
        -1/Nβ∑_,m,λ=± |g_λ()|^2 D_λ^(0)(,iω_n-iω_m)
        τ̂_3Ĝ( + ,iω_m)τ̂_3,

with self-consistent electronic Nambu Green's function Ĝ, decomposed into Pauli matrices τ̂_i, unrenomalized polaritonic Green's function D^(0), and fermionic Matsubara frequencies ω_n = (2 n+1) π /β and bosonic Matsubara frequencies ω_n = 2 n π /β, n ∈ℤ, and inverse temperature β = (k_B T)^-1. This amounts to the approximation that the bare phonon mode already contains the energy-shift renormalization due to electron-phonon coupling as the bare phonon frequency is taken from experimental data, and further renormalizations of the phonon polaritons due to electron-polariton coupling are small. The self-consistent computation of Σ̂ is initialized with a seed for the anomalous superconducting self-energy of 0.007 eV and a convergence criterion of 10^-6 eV. Further details can be found in Section D of the Supplementary Materials.


Science





§ ACKNOWLEDGMENTS

Discussions with H. Appel, S. Johnston, S. Latini, A. J. Millis, and L. Rademaker are gratefully acknowledged. M.A.S. acknowledges financial support by the DFG through the Emmy Noether programme (SE 2558/2-1). A. R. acknowledges financial support by the European Research Council (ERC-2015-AdG-694097), Grupos Consolidados (IT578-13), and  European Union's H2020  program under GA no. 676580 (NOMAD).




§ FIGURES AND TABLES















§ SUPPLEMENTARY MATERIALS


Text (Sections A to D)

References <cit.>.







 §.§ Relevant photon modes in cavity


In this work we consider a 2D material on a dielectric substrate in a nanocavity. We impose reflecting mirror boundary conditions with n⃗·B⃗=0 and n⃗×E⃗=0 for the magnetic B⃗ and electric E⃗ components of the photonic field, and n⃗=ẑ the surface normal. The size of the cavity in z direction is L_z. If the dielectric substrate has a very high dielectric constant, such as for SrTiO_3 at low temperature, it can be considered almost metallic and L_z is reduced accordingly in our effective description.

Assuming periodic boundary conditions in the x-y plane, we obtain for example for the vacuum electric field, obeying the wave equation ∇^2 E - 1/c^2∂^2 E/∂ t^2=0 with c the speed of light,

    E_x(x,y,z,t)    = E_1 exp(i k_x x) exp(i k_y y) sin(k_z z) exp(-i ω_phot(k⃗) t), 
    
    E_y(x,y,z,t)    = E_2 exp(i k_x x) exp(i k_y y) sin(k_z z) exp(-i ω_phot(k⃗) t), 
    
    E_z(x,y,z,t)    = E_3 exp(i k_x x) exp(i k_y y) cos(k_z z) exp(-i ω_phot(k⃗) t),

with ω_phot(k⃗) = c |k⃗|, and 

    k_x    = 2 π l/L_x,   l ∈ℕ_0
    
    k_y    = 2 π m/L_y,   m ∈ℕ_0
    
    k_z    = π n/L_z  n ∈ℕ_0.

We assume L_x and L_y to be large to obtain a fine momentum grid in the x-y plane. By contrast L_z is assumed to be small (L_z ≪ L_x, L_y), implying that for n = 1 the photon energy is at least c π/L_z well above typical phonon energy scales and thus irrelevant to the problem of our interest. We retain only the n=0, k_z=0 component that has constant mode amplitude along the z direction. Thus we will use only one mode for each in-plane momentum q⃗=(q_x,q_y) with

    E_x(x,y,z,t)    = 0,  
    
    E_y(x,y,z,t)    = 0, 
    
    E_z(x,y,z,t)    = E_3 exp(i q_x x) exp(i q_y y) exp(-i ω_phot(k⃗) t).





 §.§ Phonon-photon Hamiltonian


We consider the generic Hamiltonian for phonon-photon coupling <cit.>,

    H_phon-phot   = H_0 + H', 
    
    H_0    = Ω∑_q⃗ b^†_q⃗ b^_q⃗ + ∑_q⃗ω_phot(q⃗) a^†_q⃗ a^_q⃗, 
    
    H'    = -e/Mc∑_j P⃗_j ·A⃗(R⃗_j) + e^2/2Mc^2∑_j A⃗(R⃗_j) ·A⃗(R⃗_j).

Throughout we approximate the phonon dispersion relevant for FeSe/SrTiO_3 with a dispersionless Ω = 92 meV <cit.>. Here q⃗ summations are over the first Brillouin zone [-π,π)^2 in the 2D square lattice with lattice constant a=1, implying a high-frequency cutoff to the photons, which is irrelevant to the electron-boson physics happening at much lower energy. For the photon, we take only the mode polarized along the ẑ direction parallel to the phonon dipoles, and restrict it to the lowest branch q_z=0 due to cavity confinement as discussed above, implying ω_phot(q⃗) = c |q⃗| = c √(q_x^2+q_y^2).

We write the phononic dipole current operator via bosonic operators

    J⃗_j    ≡e/MP⃗_j = i e ∑_q⃗(Ω/2 N M)^1/2ξ̂_q⃗(b^†_q⃗ - b^_-q⃗) e^-i q⃗R⃗_j≡∑_q⃗1/√(N)J⃗(q⃗) e^-i q⃗R⃗_j,

with polarization vector ξ̂_q⃗ = ẑ, and similarly for the relevant z component of the photonic vector potential 

    A_z (R⃗_j)    ≡∑_q⃗(2 π c^2/ω_phot(q⃗) ν_0)^1/2(a^†_q⃗ + a^_-q⃗) e^-i q⃗R⃗_j≡∑_q⃗c/√(ν_0) A_μ(q⃗) e^-i q⃗R⃗_j,

assuming periodic boundary conditions inside the 2D plane. Here b^†_ (b^†_) creates (annihilates) a phonon with
wavevector ; a^†_ (a^†_) creates (annihilates) a cavity photon with wavevector . N is the number of unit cells, V the system volume, ν_0 ≡ V/N the unit cell volume, and e and M the ionic charge and reduced mass, respectively, related to the relative motion of positively and negatively charged ions in the optical phonon mode. In momentum space we have

    J_z(q⃗) ≡ i e (Ω/2 M)^1/2(b^†_q⃗ - b^_-q⃗), 
    
    A_z (q⃗) ≡(2 π/ω_phot(q⃗))^1/2(a^†_q⃗ + a^_-q⃗).


Now we first diagonalize the bare photon plus A^2 terms of the Hamiltonian,

    H_0,phot   = ∑_q⃗ω_phot(q⃗) a^†_q⃗ a^_q⃗
       = 1/2∑_q⃗(P_A,q⃗ P_A,-q⃗ + ω_phot(q⃗)^2 X_A,q⃗ X_A,-q⃗), 
    
    H_A^2   = 1/2∑_q⃗ω_P^2 X_A,q⃗ X_A,-q⃗,

Here we introduced canonical position and momentum operators for photon degrees of freedom,

    X_A,q⃗   ≡√(1/2 ω_phot(q⃗) )(a^_q⃗ + a^†_-q⃗), 
    
    P_A,q⃗   ≡ -i √(ω_phot(q⃗)/2)(a^_-q⃗ - a^†_q⃗).

We also defined the phononic plasma frequency

    ω_P   ≡√(4 π e^2/M ν_0) = √(4 π e^2/M ν_0, 2D L_z),

which for the 2D system in the cavity is governed by the length of the vacuum inside the cavity in z direction, L_z, and the 2D unit cell area ν_0, 2D. The expressions above are given in cgs units. In the SI system, ω^SI_P = √(e^2/M ϵ_0 ν_0, 2D L_z) with the vacuum permittivity ϵ_0.












The bilinear J · A coupling term is written as 

    H_J · A   = - 1/√(ν_0)∑_q⃗J⃗(q⃗) ·A⃗(-q⃗) 
       = -∑_q⃗ω_P X_A,q⃗ P_B,q⃗,


where it is convenient to introduce canonical position and momentum operators for the phonons,

    X_B,q⃗   ≡√(1/2 Ω)(b^_q⃗ + b^†_-q⃗), 
    
    P_B,q⃗   ≡ -i √(Ω/2)(b^_-q⃗ - b^†_q⃗).

Written in these operators, the bare phonon term H_0,phon≡Ω∑_q⃗ b^†_q⃗ b^_q⃗ takes the form

    H_0,phon   = 1/2∑_q⃗(P_B,q⃗ P_B,-q⃗ + Ω^2 X_B,q⃗ X_B,-q⃗).


The total phonon-photon Hamiltonian is now written as pairs of coupled harmonic oscillators,

    H_phon-phot   = H_0,phot + H_0,phon + H_A^2 + H_J · A
       = 1/2∑_( P_A,q⃗ P_A,-q⃗ + P_B,q⃗ P_B,-q⃗ + (ω_phot(q⃗)^2 + ω_P^2) X_A,q⃗ X_A,-q⃗ + 
                   + Ω^2 X_B,q⃗ X_B,-q⃗ - 2 ω_P X_A,q⃗ P_B,q⃗).

In order to diagonalize this Hamiltonian, we introduce a transformation

    P̃_B,   ≡Ω X_B,, 
    X̃_B,   ≡ -Ω^-1 P_B,,

which leaves the canonical commutator unchanged but interchanges position and momentum operators. 
The phonon-photon Hamiltonian is then compactly represented as

    H_phon-phot   = 1/2∑_q⃗[
    [  P_A,q⃗; P̃_B,q⃗ ]]^T
    [
    [ 1 0; 0 1 ]]
    [
    [  P_A,-q⃗; P̃_B,-q⃗ ]] + 
            + 1/2∑_q⃗[
    [  X_A,q⃗; X̃_B,q⃗ ]]^T
    [
    [ ω_phot(q⃗)^2 + ω_P^2                 Ωω_P;                 Ωω_P                  Ω^2 ]]
    [
    [  X_A,-q⃗; X̃_B,-q⃗ ]].

Diagonalization is now achieved with the following unitary transformation to polariton canonical position and momentum operators,

    [
    [ X_+,q⃗; X_-,q⃗ ]]    = [
    [  cos(θ_)  sin(θ_); -sin(θ_)  cos(θ_) ]]
    [
    [  X_A,q⃗; X̃_B,q⃗ ]],
    
    [
    [ P_+,q⃗; P_-,q⃗ ]]    = [
    [  cos(θ_)  sin(θ_); -sin(θ_)  cos(θ_) ]]
    [
    [  P_A,q⃗; P̃_B,q⃗ ]],

which leaves canonical commutation relations intact. The resulting phonon-photon Hamiltonian expressed in polaritonic operators is 

    H_phon-phot   = 1/2∑_q⃗[
    [ P_+,q⃗; P_-,q⃗ ]]^T
    [
    [ 1 0; 0 1 ]]
    [
    [ P_+,-q⃗; P_-,-q⃗ ]] + 
            + 1/2∑_q⃗[
    [ X_+,q⃗; X_-,q⃗ ]]^T
    [
    [ ω_+()^2       0;       0 ω_-()^2 ]]
    [
    [ X_+,-q⃗; X_-,-q⃗ ]],

with polaritonic dispersions ω_±() fulfilling

    ω_±()^2    = 1/2( ω_phot(q⃗)^2 + ω_P^2 + Ω^2 ±√((ω_phot(q⃗)^2 + ω_P^2 + Ω^2)^2 - 4 ω_phot(q⃗)^2 Ω^2
     )).

In particular, in the long-wavelength limit one obtains

    ω_+(q⃗→ 0)    →√(Ω^2 + ω_P^2 ), 
    ω_-(q⃗→ 0)    → 0,

as shown for the semiclassical polariton dispersions in Mahan <cit.>.
The diagonalization condition is given by 

    arctan(θ_)    = ω_phot(q⃗)^2 + ω_P^2 - Ω^2 + √((ω_phot(q⃗)^2 + ω_P^2 + Ω^2)^2 - 4 ω_phot(q⃗)^2 Ω^2 )/2 Ωω_P.

Defining bosonic operators for the upper (λ=+) and lower (λ=-) polariton branches,

    X_λ,q⃗   ≡√(1/2 ω_λ() )(α^_q⃗,λ + α^†_-q⃗,λ), 
    
    P_λ,q⃗   ≡ -i √(ω_λ()/2)(α^_-q⃗,λ - α^†_q⃗,λ).

we rewrite the phonon-photon Hamiltonian in a very compact polaritonic form:

    H_phon-phot   = ∑_,λ=±ω_λ() α^†_,λ α^_,λ.


The transformation from the initial phononic degrees of freedom to the final polaritonic ones is then given by

    X_B,q⃗   = 1/Ω (sin(θ_) P_+,q⃗ + cos(θ_) P_-,q⃗).

For the bosonic operators, this implies 

    b^†_ + b^†_-   = -i sin(θ_) √(ω_+(q⃗)/Ω) (α^_-,+ -  α^†_,+) -i  cos(θ_) √(ω_-(q⃗)/Ω) (α^_-,- -  α^†_,-)),

which will give the transformation from electron-phonon to electron-polariton coupling in the following.























 §.§ Electron-polariton Hamiltonian


The electron-polariton model Hamiltonian for FeSe/SrTiO_3 inside the cavity reads

    H    = H_e-phon + H_phon-phot, 
    
      H_e-phon   = ∑_,σϵ^†_ c^†_,σc^†_,σ
        + 1/√(N)∑_,,σ g(,)c^†_+,σ
        c^†_,σ (b^†_ + b^†_- ).

Here, c^†_,σ (c^†_,σ) creates (annihilates) an electron with
wavevector  and spin σ; ϵ_ is the electronic band dispersion measured relative to the chemical potential
μ; g(,) is the momentum dependent electron-phonon coupling. The direct electron-photon coupling of electrons in the FeSe plane to the photon branch of interest is neglected, which amounts to the assumption that the paramagnetic electronic current density j⃗ inside the FeSe layer is perfectly two-dimensional, thus not coupling to the photonic vector potential A⃗ which points perpendicular to the plane, implying j⃗·A⃗≈ 0.

Adopting the FeSe/SrTiO_3 single-band model from Rademaker et al. <cit.>, we take an electronic band dispersion ϵ_ = -2t[cos(k_x a)+cos(k_y a)] - μ, where a is the in-plane lattice constant. We set t = 0.075eV and use as an initial guess μ = -0.235eV, which is adjusted during the self-consistent calculations (see below) to a fixed band filling n_↑ = n_↓ = 0.07 for each spin. We neglect the fermion momentum dependence in the electron-phonon coupling g(,)=g(),
where  is the momentum transfer, and use g()=g_0 exp(-||/q_0). Here, g_0 is adjusted to fix the total
dimensionless coupling strength λ≈ 0.18 of the electron-phonon interaction in absence of the cavity coupling, and q_0 sets the range of the interaction in momentum
space. 

The electron-polariton expressed in polaritonic bosonic operators is obtained via Eq. (<ref>) as

    H     = ∑_,σϵ^†_ c^†_,σc^†_,σ
        + 1/√(N)∑_,,σ,λ=± c^†_+,σ
        c^†_,σ (g^*_λ() α^†_-,λ + g_λ() α^†_,λ) 
        + ∑_q⃗, λ=±ω_λ(q⃗) α^†_q⃗,λα^_q⃗,λ,

where 

    g_+()    = i sin(θ_) √(ω_+(q⃗)/Ω)  g_0 exp(-||/q_0), 
    
    g_-()    = i cos(θ_) √(ω_-(q⃗)/Ω)  g_0 exp(-||/q_0).

The couplings are thus fully determined through Eqs. (<ref>, <ref>) in connection with Eqs. (<ref>) and (<ref>).
The polariton branches and couplings to the electrons are shown in Fig. <ref> in the main text.




 §.§ Migdal-Eliashberg simulations


The electronic self-energy in Migdal-Eliashberg theory on the Matsubara frequency axis employing Nambu notation reads <cit.>

    Σ̂(,iω_n) =
        iω_n[1-Z(,iω_n)]τ̂_0 + χ(,iω_n)τ̂_3 +
        ϕ(,iω_n)τ̂_1,

where τ̂_i are the Pauli matrices, Z(,iω_n) and χ(,iω_n) renormalize
the electronic single-particle mass and band dispersion, respectively, and ϕ(,iω_n) is the anomalous
self-energy, which vanishes in the normal state. In Migdal-Eliashberg theory, the self-energy corresponding to the Hamiltonian (<ref>) is computed by
self-consistently evaluating 

    Σ̂(,iω_n) =
        -1/Nβ∑_,m |g()|^2 D^(0)(,iω_n-iω_m)
        τ̂_3Ĝ( + ,iω_m)τ̂_3,

where D^(0)(,iω_ν) = -2Ω/Ω^2 + ω_ν^2 is the bare phonon
propagator, Ĝ^-1(,iω_n) = iω_nτ̂_0 - ϵ_τ̂_3 -
Σ̂(,iω_n) is the dressed electron propagator, N is number of momentum grid points, and
β=1/(k_BT) is the inverse temperature.

Inside the cavity with ω_P>0, these well-known equations are modified to account for the Hamiltonian (<ref>) by using polariton branches λ=± instead of the phonon:

    Σ̂(,iω_n) =
        -1/Nβ∑_,m,λ=± |g_λ()|^2 D_λ^(0)(,iω_n-iω_m)
        τ̂_3Ĝ( + ,iω_m)τ̂_3,

where D_λ^(0)(,iω_ν) = -2ω_λ()/ω_λ()^2 + ω_ν^2 is the bare polariton propagator,

In practice, we use an initial guess of 0.007 eV for the anomalous self-energy and run the self-consistency until a convergence to better than 10^-6 eV is achieved. The 2D momentum grid to sample the Brillouin zone is chosen as 2000 × 2000 and convergence checked by comparing against 4000 × 4000 grids in selected cases. For the patch around q=0 we avoid the point q=0 where the lower polariton branch becomes soft since the corresponding propagator diverges in the static ω_ν=0 case. Under the q integral this divergence is cured. We therefore apply a q coarse graining by averaging 1/N_small∑̃_q |g()|^2 D^(0)(,iω_ν) over N_small small patches (∑̃_q is the sum inside the momentum patch around q=0), and using this averaged function in lieu of |g(0)|^2 D^(0)(0,iω_ν), again checking convergence in the momentum grid. The momentum convolution in Equations (<ref>) and (<ref>) is performed by fast Fourier transforms to a real-space grid and products on the real-space grid. The Matsubara cutoff is 0.4 eV for the frequency summations, and convergence in this cutoff also checked. 





