Weak and strong-coupling dynamics of quantum emitters in the vicinity of metallic films: role of electron spill-out
V. Karanikolas, I. Thanopulos, J.D. Cox, T. Kuroda, J. Inoue, N.A. Mortensen, E. Paspalakis, C. Tserkezis
aa r X i v : . [ phy s i c s . op ti c s ] F e b Weak and strong-coupling dynamics of quantum emitters in the vicinity of metallicfilms: role of electron spill-out
Vasilios Karanikolas a , ∗ Ioannis Thanopulos b , Joel D. Cox c,d , Takashi Kuroda e ,Jun-ichi Inoue e , N. Asger Mortensen c,d , Emmanuel Paspalakis b , and Christos Tserkezis c † a International Center for Young Scientists (ICYS),National Institute for Materials Science (NIMS) 1-1 Namiki,Tsukuba, Ibaraki 305-0044, Japan, b Materials Science Department,School of Natural Sciences, University of Patras,Patras 265 04, Greece, c Center for Nano Optics,University of Southern Denmark, Campusvej 55, DK-5230 Odense M,Denmark, d Danish Institute for Advanced Study,University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark, e National Institute for Materials Science (NIMS) 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan (Dated: February 23, 2021)The relaxation of a quantum emitter (QE) near metal-dielectric layered nanostructures is inves-tigated, with focus on the influence of plasmonic quantum effects. The Green’s tensor approach,combined with the Feibelman d -parameter formalism, is used to calculate the Purcell factor and thedynamics of a two-level QE in the presence of the nanostructure. Focusing on the case of Na, weidentify electron spill-out as the dominant source of quantum effects in jellium-like metals. Our re-sults reveal a clear splitting in the emission spectrum of the emitter, and non-Markovian relaxationdynamics, implying strong light–matter coupling between them, a coupling that is not prevented bythe quantum-informed optical response of the metal. Introduction — Confining light in sub-diffraction vol-umes via surface plasmon polaritons (SPPs) has becomecommon practice to enhance light–matter interactionsin the past decade [1]. In information-processing ap-plications, plasmonic nanostructures are combined withmulti-level quantum emitters (QEs) which, dependingon the desired functionality, can be natural (atoms,molecules) or artificial (quantum wells and dots, de-fects in diamonds, collective states in transition metaldichalcogenides or hexagonal boron nitride) [2–7]. Byanalogy to cavity QED [8, 9], metallic nanostructuresplay in such cases the role of the cavity [10–13]. Ad-vanced fabrication techniques now enable the engineer-ing of metal cavities with characteristic dimensions be-low nm, and the precise positioning of QEs withinthem, to produce promising templates for bright single-photon emitters [14–18]. The plasmonic cavity is typ-ically formed by a metallic substrate interacting withanother metallic film or with dropcasted nanoparticles;separating the two plasmonic components by a thin di-electric spacer containing QEs in the form of defects,molecules, or quantum dots [19–22], the resulting intenseelectromagnetic (EM) fields, confined within minimizedvolumes, trigger extreme light–matter interactions.Traditionally, in the local-response approximation(LRA) [23], the optical response of metals is describedthrough the Drude model, or through experimentallymeasured bulk permittivities. But, as plasmonic cavi-ties become narrower, the interaction strengths with thetherein confined QEs are overestimated by LRA, as com-pared to experimental results [24–26], calling for amend-ment of the theoretical models [27–31]; the main miss-ing element is information about quantum effects in the metal. As a first extension of LRA, the hydrodynamicDrude model (HDM), which describes the motion of thecompressible free-electron gas as a convective fluid, hasmet with considerable success [32–35], yet effects associ-ated with electron spill-out still require a self-consistenttreatment [36–38]. However, the ultimate goal of an ab initio optical response calculation for structures withsizes of nm order at room temperature is out ofcomputational reach [39]. A more efficient approach isbased on the work of Feibelman [40], who bridged EMand ab initio calculations through appropriate surface-response functions, the d parameters, which need be cal-culated only once for a given surface, and can be im-plemented in a quantum-informed surface-response for-malism (SRF) [31, 41–43]. The full QED problem isthen described by quantum-corrected mesoscopic bound-ary conditions at the metal surface [31, 41]. MacroscopicQED [44, 45] allows to directly apply the rigorous Green’sfunction approach developed within LRA, now correctedby SRF, without invoking quantum effects in a quasi-normal mode formalism [46–48].To date, focus has been cast on the interaction ofQEs with plasmonic nanostructures described by non-local dielectric functions (e.g. within HDM), identifyingsignificant differences in both weak [49, 50] and strong-coupling regimes [51, 52], as compared to LRA. Here,we go one step further to explore, for the first time, therole of quantum-informed plasmonics in the populationdynamics of the QE excited state under strong-couplingconditions, where QE and environment coherently ex-change energy (observable through Rabi oscillations inthe population and Rabi energy splittings in the emis-sion spectra). We focus on insulator/metal (IM), insu-lator/metal/insulator (IMI), and metal/insulator/metal(MIM) geometries [33, 53], for which the strong-couplingregime proves reachable even when quantum correctionsin the metal are taken into account. Theory — The QE is approximated as a two-level sys-tem with ground state | g i and excited state | e i , transi-tion frequency ω , and dipole moment µ . Initially, theQE is in the excited state, and the EM field is in its vac-uum state | i i = | e i ⊗ | i , while the final state is | f i = | g i⊗ ˆ f † i ( r , ω ) | i , where the QE relaxes to the ground stateby emitting a photon or exciting SPPs dressed statessupported by the metallic nanostructures [54], and ˆ f † de-notes the bosonic creation operator of the dressed state i .The QE relaxation rate Γ is given by Γ( r , ω ) = Γ ˜Γ( r , ω ) ,where Γ = ω µ / πc ~ ε is the vacuum rate and ˜Γ( r , ω ) is the Purcell factor of a QE placed at r = (0 , , z ) abovea metal/dielectric interface; for a transition dipole mo-ment perpendicular to the interface, the Purcell factorhas the form ˜Γ z ( r , ω ) = √ ε d − c ω Re ( ∞ Z d k s k s R TM k d,z k d e ik d,z z ) , (1)where R TM is the generalized transverse magnetic (TM)Fresnel coefficient [55, 56].Our material of choice is Na, whose work function islower than that of the common plasmonic metals, Au andAg, so as to focus on the electron spill-out effects [57,58]. For the single dielectric/metal interface, the TMreflection coefficients take the form (assuming d k = 0 ,valid for Na) [43] R TM = ε m k d ,z − ε d k m ,z + i ( ε m − ε d ) k s d ⊥ ( ω ) ε m k d ,z + ε d k m ,z − i ( ε m − ε d ) k s d ⊥ ( ω ) , (2)where ε d and ε m denote the permittivity of the dielectric(d) and metal (m), respectively, k = ω/c is the free-space wave vector, k j = √ ε j k (with j =m,d) is the wavevector of each medium, analyzed in in-plane ( k s ) andnormal ( k j,z = q k j − k s ) components. For Na we usea Drude model ε ( ω ) = 1 − ω / (cid:0) ω + iωγ (cid:1) , with ω p be-ing the plasma frequency and γ the damping rate, takenas ~ ω p = 5 . eV and ~ γ = 0 . eV, while d ⊥ is given byLorentzian-fitted data extracted from ab initio calcula-tions for jellium [42, 43].The relaxation of the quantum emitter is described byits emission spectrum [59, 60] S ( ω, r ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( µ ω / √ πε c ) ˆn · G ( ω, r , r d ) ω − ω − R ∞ dω ′ J ( ω , ω ′ , r ) / ( ω ′ − ω ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (3)in the frequency domain, and by the time-dependent QEexcited state probability amplitude, c ( t ) , given by thesolution of the integro-differential equation [62–65] dc ( t ) dt = i Z t dt ′ K ( t − t ′ ) c ( t ′ ) . (4) Figure 1. A QE with z -oriented transition dipole moment,placed nm away from an air/Na IM interface. (a) Purcellfactor and (b) emission spectrum within LRA (black dashedlines) and SRF (red solid lines). The kernel of Eq. (4) is given by K ( τ ) = ie iω τ R ∞ J ( ω , ω, r ) e − iωτ dω , where J ( ω , ω, r ) =˜Γ i ( ω, r ) / [2 πτ ] is the spectral density, ω is the energydifference between the ground and excited QE states,and τ = 1 / Γ ( ω ) is its free-space lifetime [62–65]. InEq. (3), r d is the position where the signal is detected,and r is the position of the QE, while ω is the emissionfrequency of the combined QE/Na system [59, 60].Eq. (3) implies that if the QE–nanostructure couplingstrength is high enough, the emission spectrum willfeature a doublet of emission peaks, whose energydifference defines the Rabi splitting ~ Ω . More details forthe macroscopic QED model used can be found on [61]. Results — We begin our investigation of the influenceof spill-out on the emission properties of QEs by revisit-ing the simple IM case. In Fig. 1(a) we present the Pur-cell factor as a function of the QE transition energy fora QE placed nm away from an air–Na IM interface, asevaluated by LRA and SRF. We directly observe that thehighest Purcell factor value within SRF has a redshiftedpeak value that is one order of magnitude smaller thanthe corresponding LRA result, and the SPP-originatedresonance around ~ ω = ~ ω p / √ . eV is broader.These features are related to additional material lossesdue to surface-enabled Landau damping [66]. Neverthe-less, at lower frequencies, away from the SPP resonance,the SRF Purcell factor is in fact higher than the onepredicted by LRA. Similar results have been discussedin Ref. 43, thus validating the Green’s tensor formalismemployed here.The emission spectrum of a QE with transition energy ~ ω = 4 . eV and free-space lifetime τ = 0 . ns is pre-sented in Fig. 1(b), where the Purcell factors of Fig. 1(a)are used. In the LRA description, the emission spec-trum features two sharp peaks, splitting by ~ Ω = 1 . eVaround the QE transition energy. When the quantum as-pects of the response of the metal are introduced throughSRF, the two emission peaks become broader due to thehigher material losses, but the Rabi splitting increasesto ~ Ω = 1 . eV; the latter phenomenon stems from thedependence of the QE emission on the entire spectrum of Figure 2. Purcell factor within (a) LRA and (b) SRF of a QEplaced nm away from an IMI nanostructure, for differentmetal thicknesses D = 10 (black solid line), 20 (red dottedline) and nm (blue dashed line). The thin purple showsthe corresponding IM result. the Purcell factor, which is higher within SRF away fromthe SPP resonance [Fig. 1(a)] [42].Turning to the IMI case, where air is considered asthe insulating medium, we examine the correspondingPurcell factor in Figs. 2(a) and (b) within LRA and SRF,respectively, for metal layer thicknesses D = 10 , , and nm. Specifically, we set the distance between the QEand the insulator/metal interface to z QE = 5 nm, notingthat the Purcell factors calculated for the IM and theIMI geometry exhibit negligible differences for smallervalues z QE . nm. In the case of LRA, the Purcellfactor has its highest value at the SPP energy for all layerthicknesses. On the other hand, when SRF is taken intoaccount, the maximum Purcell factor is roughly halved,and its peak positions redshift significantly as the metalthickness decreases. For the thinnest layer ( D = 10 nm),the resonance within SRF has shifted by more than ~ ω =0 . eV compared to LRA. For increasing metal thickness(above D = 50 nm), the Purcell factor for the IMI caseconverges to the IM result.While the IMI architecture is instructive, we will fo-cus on the MIM structure, which is expected to displayall the advantages of plasmonic cavities discussed in theintroduction. In Figs. 3(a,b) we compare LRA and SRFcontour plots of the Purcell factor as a function of theQE emission energy and the thickness D of the insulatorlayer (assumed air here) for a QE placed in the middleof the insulator layer of the MIM structure ( z QE = 0 ).In both plots we observe that, for small D , the Purcellfactor is significantly enhanced compared to the referencefree-space value. In LRA, the peak of the Purcell factorcoincides with the SPP resonance, while for higher ener-gies there is a sharp drop of its value, connected with theabsence of plasmon polariton modes in the band-gap re-gion of the MIM structure [67]. On the other hand, SRFreproduces the expected redshift of the Purcell factorspectrum, while considerable enhancement is observedinside the band-gap, as a result of the quantum effectscaptured by the Feibelman parameters. As the thicknessof the insulator increases, surface effects become less im- Figure 3. Contour plot of the Purcell factor of a QE, as afunction of emission energy and dielectric layer thickness D ,calculated within (a) LRA and (b) SRF, and the correspond-ing contours of the emission spectra as a function of the QEemission and transition energies obtained with (c) LRA and(d) SRF. The QE, with τ = 0 . ns, is placed in the middle ofthe insulator layer of an MIM structure with dielectric-layerthickness D = 4 nm.Figure 4. Population of the excited state of a two-level QEwith ~ ω = 3 eV and τ = 0 . ns, as a function of time t . (a)The QE is placed nm away from an IM geometry. (b) TheQE is placed in the middle of a MIM geometry with dielectric-layer thickness D = 2 nm. In both panels, black and red linescorrespond to results within LRA and SRF, respectively. portant and the spectra obtained within the two modelsconverge. What is evident from the large Purcell-factorenhancements of Figs. 3(a,b) is that the QE–MIM inter-action can enter the strong coupling regime, as verified bythe emission-spectrum contour plots of Fig. 3(c,d), whereclear anticrossings can be observed not only within LRA,but also in the SRF case. Naturally, the emission double-peak features are sharper in the LRA description, and arebroadened when additional loss channels due to surfaceeffects are considered, but the observed Rabi splittingclearly survives.The QE excited state population dynamics presentedin Fig. 4 exhibits strong oscillations, with very differentfeatures from simple exponential relaxation within theMarkovian (weak coupling) approximation, where mem-ory effects in Eq. (4) are ignored. For a QE with tran- Figure 5. Population of the excited state of a two-level QEas a function of time, within LRA (a,c) and SRF (b,d). TheQE is placed in the middle of a MIM geometry of thickness D = 2 nm. (a,b) The values of τ considered are displayed inthe insets. The transition energy of the QE is ~ ω = 4 . eV(LRA) and ~ ω = 3 . eV (SRF). (c,d) The values of ~ ω considered are shown in the insets, for a QE with τ = 0 . ns. sition energy ~ ω = 3 eV and vacuum relaxation rate τ = 0 . ns, Fig. 4(a) shows the excited state dynam-ics when z QE = 1 nm in the IM geometry. In the LRAdescription, we observe that after a few rapid oscillations,the QE relaxes to the ground state within a time-span of fs. Faster relaxation to the ground state is predictedin SRF, although a closer look in the inset of Fig. 4(a) re-veals that in this case the excited state remains partiallypopulated, | c nL | = 0 . , even at longer times.To directly compare IM and MIM geometries, we con-sider a QE centered in an insulating (air) layer of D =2 nm thickness sandwiched between metal regions, i.e.,maintaining a nm distance from the metal/dielectricinterfaces. Here, the QE excited state population dy-namics presented in Fig. 4(b) is characterized by strongoscillatory behavior, with large population values persist-ing over long time spans. While the QE eventually fullyrelaxes, at later times, to the ground state in LRA, theSRF calculations accounting for the quantum responsepredicts that about ∼ % of its initial population re-mains trapped in the excited state.In Fig. 5 we present the excited state population dy-namics of a QE placed in the middle of a MIM cavityhaving an insulator (air) layer thickness of D = 2 nmwhen different free-space relaxation rates τ are con-sidered. Choosing the associated QE transition energymatching the highest value shown in the Purcell factorspectra of Fig. 3, LRA and SRF models are contrastedin Figs. 5(a,b) for ~ ω = 4 . eV and ~ ω = 3 . eV, re-spectively. In both panels we observe strong and rapid Rabi oscillations; as the value of the free-space lifetimedecreases from τ = 10 to . ns, the QE/MIM cou-pling increases, the oscillation period decreases, and non-Markovian effects become stronger, as anticipated. Wealso observe that the population oscillations are denserin the SRF case.Finally, we explore the influence of the QE transi-tion energy, ~ ω , on the excited state population dy-namics in Figs. 5(c,d), taking a QE free-space lifetime of τ = 0 . ns. Different values of the QE transition energyare considered, including those matching the SPP ener-gies indicated in the Purcell factor spectrum of Fig. 3.Fig. 5(c) shows the LRA calculation, where we observethat, after a few oscillations, the population dynam-ics remains partially trapped—although, eventually, itfully relaxes at later times. As the QE transition en-ergy decreases, moving out of resonance with the SPP atits associated Purcell enhancement, the partially tran-sient trapped population steadily decreases. In the SRFcase, we also observe a highly oscillatory behavior of theexcited state population density, a sign that the non-Markovian signature is retained. Although the partiallyexcited population trapping is smaller, strong oscillationsare still present, with reduced period, a sign that theoverall QE–MIM interaction is higher when the quantumaspects of the metal response are taken into account. Summary — We theoretically explored the emissionproperties of a QE placed in proximity to IM, IMI, andMIM geometries in the weak and strong coupling regimes.Focusing on Na as a plasmonic metal dominated by elec-tron spill-out, we described its quantum-informed opti-cal response through SRF, in which the Feibelman d pa-rameters for the centroid of induced charge incorporatefirst-principles calculations. Although quantum effectsintroduce additional losses in the metal, the QE/MIMcavity system is found to operate in the strong-couplingregime, manifested as a Rabi splitting in the emissionspectrum, and through Rabi oscillations in the relax-ation dynamics of the QE excited state, with a highlynon-Markovian character. SRF bridges, therefore, ab ini-tio approaches with classical EM calculations in a simpleand elegant way, suitable for quantum-optical studies atthe nanoscale. Acknowledgments — J. D. C. is a Sapere Aude researchleader supported by Independent Research Fund Den-mark (grant no. 0165-00051B). N. A. M. is a VILLUMInvestigator supported by VILLUM FONDEN (grant no.16498). E.P.’s work is co-financed by Greece and the Eu-ropean Union (project code name POLISIMULATOR). ∗ [email protected] † [email protected][1] E. Ozbay, Science , 189 (2006). [2] A. V. Akimov, A. Mukherjee, C. L. Yu, D. E. Chang,A. S. Zibrov, P. R. Hemmer, H. Park, and M. D. Lukin,Nature , 402 (2007).[3] D. E. Chang, A. S. Sørensen, E. A. Demler, and M. D.Lukin, Nature Physics , 807 (2007).[4] A. Goban, C.-L. Hung, S.-P. Yu, J. Hood, J. Muniz, J.Lee, M. Martin, A. McClung, K. Choi, D. Chang, O.Painter, and H. Kimble, Nature Communications , 3808(2014).[5] J. D. Hood, A. Goban, A. Asenjo-Garcia, M. Lu, S.-P.Yu, D. E. Chang, and H. J. Kimble, Proceedings of theNational Academy of Sciences , 10507 (2016).[6] A. I. Fernández-Domínguez, S. I. Bozhevolnyi, and N. A.Mortensen, ACS Photonics , 3447 (2018).[7] M. Koperski, D. Vaclavkova, K. Watanabe, T. Taniguchi,K. S. Novoselov, and M. Potemski, Proceedings of theNational Academy of Sciences , 13214 (2020).[8] L. Mandel and E. Wolf, Optical Coherence and QuantumOptics (Cambridge University Press, 1995).[9] M. O. Scully and M. S. Zubairy,
Quantum Optics (Cam-bridge University Press, 1997).[10] P. Törmä and W. L. Barnes, Reports on Progress inPhysics , 013901 (2015).[11] M. Pelton, Nature Photonics , 427 (2015).[12] S. I. Bozhevolnyi and J. B. Khurgin, Nature Photonics , 398 (2017).[13] J. J. Baumberg, J. Aizpurua, M. H. Mikkelsen, and D.R. Smith, Nature Materials , 668-678 (2019).[14] C. Manolatou and F. Rana, IEEE Journal of QuantumElectronics , 435 (2008).[15] G. M. Akselrod, C. Argyropoulos, T. B. Hoang, C. Ciracì,C. Fang, J. Huang, D. R. Smith, and M. H. Mikkelsen,Nature Photonics , 835 (2014).[16] A. Rose, T. B. Hoang, F. McGuire, J. J. Mock, C. Ciracì,D. R. Smith, and M. H. Mikkelsen, Nano Letters , 4797(2014).[17] R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A.Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess,and J. J. Baumberg, Nature ,127 (2016).[18] M. A. May, D. Fialkow, T. Wu, K.-D. Park, H. Leng,J. A. Kropp, T. Gougousi, P. Lalanne, M. Pelton, andM. B. Raschke, Advanced Quad. Tech. , 1900087 (2020).[19] S. I. Bogdanov, M. Y. Shalaginov, A. S. Lagutchev, C.-C. Chiang, D. Shah, A. S. Baburin, I. A. Ryzhikov, I. A.Rodionov, A. V. Kildishev, A. Boltasseva, and V. M.Shalaev, Nano Letters , 4837 (2018).[20] J. M. Katzen, C. Tserkezis, Q. Cai, L. H. Li, J. M. Kim,G. Lee, G.-R. Yi, W. R. Hendren, E. J. G. Santos, R. M.Bowman, and F. Huang, ACS Applied Materials & In-terfaces , 19866 (2020).[21] O. S. Ojambati, W. M. Deacon, R. Chikkaraddy, C.Readman, Q. Lin, Z. Koczor-Benda, E. Rosta, O. A.Scherman, and J. J. Baumberg, ACS Photonics , 2337(2020).[22] F. Ding, Y. Yang, R. A. Deshpande, and S. I. Bozhevol-nyi, Nanophotonics , 1129 (2018).[23] U. Hohenester, Nano and Quantum Optics: An Introduc-tion to Basic Principles and Theory (Springer, 2020).[24] C. Ciracì, R. T. Hill, J. J. Mock, Y. Urzhumov, A. I.Fernández-Domínguez, S. A. Maier, J. B. Pendry, A.Chilkoti, and D. R. Smith, Science , 1072 (2012).[25] K. J. Savage, M. M. Hawkeye, R. Esteban, A. G. Borisov,J. Aizpurua, and J. J. Baumberg, Nature , 574(2012). [26] J. A. Scholl, A. L. Koh, and J. A. Dionne, Nature ,421 (2012).[27] I. Romero, J. Aizpurua, G. W. Bryant, F. J. García deAbajo, Optics Express , 9988 (2006).[28] R. Esteban, A. G. Borisov, P. Nordlander and J. Aizpu-rua, Nature Communications , 825 (2012).[29] Y. Luo, A. I. Fernández-Domínguez, A. Wiener, S. A.Maier and J. B. Pendry, Physical Review Letters ,093901 (2013).[30] S. Raza, N. Stenger, S. Kadkhodazadeh, S. V. Fischer, N.Kostesha, A.-P. Jauho, A. Burrows, M. Wubs, and N. A.Mortensen, Nanophotonics , 131 (2013).[31] Y. Yang, D. Zhu, W. Yan, A. Agarwal, M. Zheng,J. D. Joannopoulos, P. Lalanne, T. Christensen, K. K.Berggren, and M. Soljačić, Nature , 248 (2019).[32] A. Moreau, C. Ciracì, and D. R. Smith, Physical ReviewB , 045401 (2013).[33] S. Raza, T. Christensen, M. Wubs, S. I. Bozhevolnyi, andN. A. Mortensen, Physical Review B , 115401 (2013).[34] N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard,and S. I. Bozhevolnyi, Nature Communications , 3809(2014).[35] C. Tserkezis, A. T. M. Yeşilyurt, J.-S. Huang, and N. A.Mortensen, ACS Photonics , 5017 (2018).[36] G. Toscano, J. Straubel, A. Kwiatkowski, C. Rockstuhl,F. Evers, H. Xu, N. A. Mortensen and M. Wubs, NatureCommunications , 7132 (2015).[37] W. Yan, Physical Review B , 115416 (2015).[38] C. Ciracì and F. Della Sala, Physical Review B ,205405 (2016).[39] A. Varas, P. García-González, J. Feist, F. García-Vidal,and A. Rubio, Nanophotonics , 409 (2016).[40] P. Feibelman, Progress in Surface Science , 287 (1982).[41] W. Yan, M. Wubs, and N. Asger Mortensen, PhysicalReview Letters , 137403 (2015).[42] T. Christensen, W. Yan, A.-P. Jauho, M. Soljačić, andN. A. Mortensen, Physical Review Letters , 157402(2017).[43] P. A. D. Gonçalves, T. Christensen, N. Rivera, A.-P.Jauho, N. A. Mortensen, and M. Soljačić, Nature Com-munications , 366 (2020).[44] N. Rivera and I. Kaminer, Nature Reviews Physics , 538(2020).[45] K. Head-Marsden, J. Flick, C. J. Ciccarino,and P. Narang, Chemical Reviews (2021),10.1021/acs.chemrev.0c00620.[46] P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P.Hugonin, Laser & Photonics Reviews , 1700113 (2018).[47] M. K. Dezfouli, C. Tserkezis, N. A. Mortensen, and S.Hughes, Optica , 1503 (2017).[48] S. Franke, S. Hughes, M. K. Dezfouli, P. T. Kristensen,K. Busch, A. Knorr, and M. Richter, Physical ReviewLetters , 213901 (2019).[49] C. Tserkezis, N. Stefanou, M. Wubs, and N. A.Mortensen, Nanoscale , 17532 (2016).[50] T. Christensen, W. Yan, S. Raza, A.-P. Jauho, N. A.Mortensen, and M. Wubs, ACS Nano , 1745 (2014).[51] C. Tserkezis, M. Wubs, and N. A. Mortensen, ACS Pho-tonics , 133 (2018).[52] C. Tserkezis, A. I. Fernández-Domínguez, P. A. D.Gonçalves, F. Todisco, J. D. Cox, K. Busch, N. Stenger,S. I. Bozhevolnyi, N. A. Mortensen, and C. Wolff, Re-ports on Progress in Physics 83, 082401 (2020).[53] C. L. C. Smith, N. Stenger, A. Kristensen, N. A. Mortensen, and S. I. Bozhevolnyi, Nanoscale , 9355(2015).[54] H. T. Dung, L. Knöll, and D.-G. Welsch, Physical ReviewA , 043813 (2002).[55] C. T. Tai, Dyadic Green Functions in ElectromagneticTheory (Oxford University Press, Oxford, 1994).[56] S. Scheel and S. Y. Buhmann, Acta Physical Slovaca
675 (2009).[57] F. J. García de Abajo, Journal of Physical Chemistry C , 17983 (2008).[58] A. R. Echarri, P. A. D. Gonçalves, C. Tserkezis, F. J.García de Abajo, N. A. Mortensen, and J. D. Cox,arXiv:2009.10821.[59] C. Van Vlack, P. T. Kristensen, and S. Hughes, PhysicalReview B , 075303 (2012).[60] V. Karanikolas, I. Thanopulos, and E. Paspalakis, Phys-ical Review Research , 033141 (2020). [61] See Supplemental Material for details on the calculationof the Green’s tensor for the IM, IMI and MIM geome-tries including the SRF boundury conditions. The macro-scopic QED theory is used to extract Eqs. (3) and (4).[62] A. Gonzalez-Tudela, P. A. Huidobro, L. Martín-Moreno,C. Tejedor, and F. J. García-Vidal, Physical Review B , 041402(R) (2014).[63] R.-Q. Li, D. Hernángomez-Peréz, F. J. García-Vidal,and A. I. Fernández-Domínguez, Physical Review Let-ters , 107401 (2016).[64] I. Thanopulos, V. Yannopapas, and E. Paspalakis, Phys-ical Review B , 075412 (2017).[65] I. Thanopulos, V. Karanikolas, N. Iliopoulos, and E. Pas-palakis, Physical Review B , 195412 (2019).[66] C. Tserkezis, N. A. Mortensen, and M. Wubs, PhysicalReview B , 085413 (2017).[67] C. A. Marocico and J. Knoester, Physical Review A84