Quantum Many-Body Effects in the Optoelectronic Response of Plasmonic Nanostructures and their Coupling to Quantum Emitters

EUSKAL HERRIKO UNIBERTSITATEA THE UNIVERSITY OF THE BASQUE COUNTRY Department of Electricity and Electronics CAMPUS OF INTERNATIONAL EXCELLENCE Quantum Many-Body Effects in the Optoelectronic Response of Plasmonic Nanostructures and their Coupling to Quantum Emitters Thesis by Antton Babaze Aizpurua Supervised by Prof. Javier Aizpurua Iriazabal and Dr. Rubén Esteban Llorente Donostia-San Sebastián, September 2022 (cc)2022 ANTTON BABAZE AIZPURUA (cc by-sa 4.0) Ama eta Aitari, tesi honen oinarri sendoak ezartzeagatik. ACKNOWLEDGMENTS Han pasado casi cinco años desde el 2 de octubre de 2017, primer día que acudía a la oficina de Rubén, por aquel entonces en el DIPC. Javi estaba en algún viaje de los suyos. Aquel día empezaba mi primer trabajo “serio”, donde me pagaban por “estudiar”. Me pagaban poco para mi gusto, pero eso es otro tema. Ese primer día Rubén me imprimió unos 10 artículos, no sé cuántas hojas en total (muchas), “para que vayas leyendo y haciéndote una idea”, me decía. Antes de imprimirlos, claro, me explicó lo mejor que pudo de qué iban esos artículos y cómo se relacionaban con lo que iba a ser mi tesis. Tesis que culmino ahora, cinco años más tarde, con estas palabras de agradecimiento. También me mostró uno por uno los posters que había colgados en la pared del pasillo en el CFM. Creo que fingí bastante bien que más o menos estaba siguiendo lo que me decía, pero he de confesar que, en realidad, no estaba entendiendo nada. No sé si Ruben era consciente de ello, puede que sí. Tardé unos cuantos meses en saber, por ejemplo, qué era un plasmón. Creo que ahora puedo afirmar que entiendo los artículos que me imprimió hace cinco años como referencia. En gran parte, gracias al constante aprendizaje que he recibido por parte de Rubén. No sabes cuánto tiempo tardé en averiguar cómo se relacionaban las pelotitas de metal iluminadas con láser con el par de muelles acoplados que oscilan siguiendo las ecuaciones que aprendimos en segundo de carrera. Además, aparte de haber sido capaz de entender aquellos artículos, me siento bastante orgulloso de haber contribuido, aunque sea mínimamente, en el campo de la plasmónica. Contribuciones que no serían posibles sin tu conocimiento, trabajo incansable, y esfuerzo por enseñarme. Gracias, Rubén. A las pocas semanas de empezar mi doctorado entró Javi en escena. Fue en un café en el DIPC. De primeras supe que iba a ser muy sencillo tratar con él como supervisor. Así fue. Por supuesto, no sabía ni que era el director del CFM ni que era un tío tan conocido en el campo de la Nanofotónica, con un índice h de no sé cuánto. Me hacía mucha gracia cuando, en las conferencias, todo el mundo me decía: “Of course I know Javier Aizpurua”. Al principio no sabía muy bien por qué. Tampoco sabía la importancia que Javi iba a tener en mi tesis, que desde luego ha sido mucha. Urte hauetan zehar sarritan izan dut buruan Pernando Amezketarraren marrazki bizidunetan behinola entzundako esaldia: “Galdera zailik ez dago, erantzunak dira zailak”. Denon artean, nire tesian hain errazak ez ziren erantzun gutxi batzuk ematea lortu dugu. Ez dago gaizki. Horretan, pisu handia izan duzu zuk, Javi. Eskerrik asko urte hauetan v irakatsitako guztiarengatik: zientzian eta, esango nuke, batez ere zientziatik kanpo irakatsitakoagatik. Punta-puntako talde batean ikertzeko aukera emateagatik. Nire ikerketa-lana maisuki gidatzeagatik. Arazo txikien aurrean beti irtenbiderik hoberena bilatzeagatik. Nire lana aurrera ateratzeko egindako esfortzuarengatik. Nire alde egindako apustuarengatik. Eskerrik asko, Javi. In May 2018, just seven months after I started with my PhD thesis, I had the opportunity to meet Andrei Borisov in ISMO, Orsay. Andrei has been my third supervisor, along with Ruben and Javi. He has been the key person in this thesis. The one who has taught me everything I know about TDDFT and quantum plasmonics. All the research in this thesis has been carried out in extremely close collaboration with Andrei. Thank you for your patience, guidance, and willingness to teach me. Willingness to discuss with me. So many discussions. So many calls. I have checked that we have exchanged 321 e-mails during the last four years. Thank you for sharing with me your WPP codes. This thesis is also yours. Only one remark: we shouldn´t refer to sodium as a simple metal. Thank you, Andrei. Si he tenido suerte con los supervisores de mi tesis, he tenido aún mucha más suerte con los compañeros de trabajo que he tenido en el CFM. Compañeros de trabajo que son ahora amigos. No recuerdo en qué tesis leí que los agradecimientos son la parte más difícil de escribir. Para nada. Me resulta verdaderamente fácil escribir estas líneas de forma honesta. Habéis sido la hostia. No creo que nunca tenga ningún trabajo donde vaya todos los días con tantas ganas a trabajar. Ha sido súper divertido haber compartido tantos momentos con vosotros. Dentro del trabajo, y por supuesto, fuera. Con el riesgo de olvidarme de alguien importante, procedo a nombrar a varios de vosotros. Empiezo por Álvaro, Alvariño, mi masajista particular. Compañero de doctorado desde el principio hasta el final. Compañero en las clases magistrales de Rubén. Cómo me aliviaba saber que tú tampoco te enterabas de nada. O de muy poco. No sé en qué punto empezamos a investigar cosas tan distintas, cuando al principio parecía que estudiábamos exactamente lo mismo: el dipolo junto a la pelotita. Ha sido un placer haber hecho todo el doctorado contigo, haber charlado sobre ciencia, y, sobre todo, haber charlado sobre tantas otras cosas fuera de la ciencia. Muchos buenos momentos fuera de la oficina. Viajes. También dentro de la oficina. Muchas risas. Muchas estupideces. Mucho tiempo “perdido”, pero con mucho gusto. Suerte con tu tesis, tío. En nada lo tienes. Continúo con los demás muchachos de la oficina. Carlos, Charly, Carlitos. Teoremán. Qué capacidad, macho. Siempre que he tenido alguna duda sobre física o matemática, ahí estabas tú con tus teoremas. Sé que muchos eran inventados, aunque nunca lo reconozcas. A ver cuándo me pones de autor en uno de tus papers. Una máquina con inkscape y blender. Cuando nadie te ve, claro. Nunca olvidaré todos esos días de charleta en la oficina cuando todo el mundo ya se había marchado a casa. O los skypes durante la pandemia. O aquellos días previos a Navidad, cuando nos quedábamos los dos solos trabajando y tú ponías villancicos a todo volumen fingiendo ser feliz. Yo me lo pasaba bien. Me alegro de que mexicanos como tú crucen el charco para quitarnos el trabajo, wey. Jon, Joneto, qué tío. Siempre dispuesto a ayudar en lo que sea, pero también a sacarle punta vi a todo lo que yo digo. Berdin du nik zer esan, beti egongo zarela kontran. Un titán. Dibertigarria izan da jokoa, benetan. Jarraitu horrela, motel. Bruno, el más listo del CFM. Seguramente también el más rápido con el coche. Haciendo TDDFT como yo, pero bien. Zorte on doktoretzarekin. Adrián, el más gracioso, otro que se quedaba haciendo guardia hasta bien tarde en la oficina del CFM. Jonathan, con ese humor tan particular, siempre “hueveando”. Mikel, Mikelats, nire Hondarribiko bizilaguna. Tipo bikaina. Denetan ona. Plazera izan da zurekin denbora guzti hau partekatzea. Zorte on zuri ere tesiarekin, aurki bukatuko duzula. Jorge, otro grande. A ver cuándo te sacas el título de euskera. Un placer los vaciles mutuos que hemos tenido, y espero seguir teniendo. Roberto, Rober, ídolo. Me he reído mucho contigo, eta asko eskertzen dut euskararekin egiten duzun ahalegina. Ez nuen espero Asturiar batekin euskaraz hitz egiten bukatuko nuenik. Iker, mi Bilbaino favorito. Gracias por haberme dejado vacilarte tanto. Lástima la última copa. Gracias a otros muchos también: Martín, Miriam, Mattin, Fernando, Josu, Joseba, Mario, Raulillo, Txemikel... Todos vosotros habéis conseguido que sea un verdadero placer venir todos los días al CFM. Lanean egindako lagunez gain, lanetik kanpo ditudan lagunak ere aipatu nahi nituzke. Lerro gutxietan bada ere. Irungo betiko lagunak. Artalekuko Hintxak. Muerte a los Geólogos. Simon-eneako Cortijeroak. Hulk. Kemenkideak. Dultzaineroak. “Noiz bukatzen duzu proiektua?”; “Zer moduz daramazu mastera?”; “Baina, tira, zurea ez da “lana-lana”, ez?”. Eskerrik asko denoi. Bukatzeko, lerro batzuk eskaini nahi dizkiet tesi honen benetako euskarri izan zaretenei: Etxekoei. Tesi hau zuena da. Zuek gabe ezingo nuke tesi hau inoiz aurrera atera. Ama, horrenbeste gauza irakasteagatik. Nik hartutako erabakiak beti babesteagatik. Nire helburuak betetzeko beti zure esku dagoen guztia egiteagatik. Ez dut sekula ahaztuko momentu zailetan egindako esfortzua. Ondo dakizu zertaz ari naizen, noski. Aita, niregan beti jarrera kritikoa bultzatzeagatik. Gauzak esfortzuarekin ateratzen direla irakasteagatik. Bizitzan bide onetik noala behin eta berriro azpimarratzeagatik. Niregan erakutsitako konfiantzagatik. Iñaki, zuregandik jasotako animo eta laudorioengatik. Nire tesiaz, eta bizitzaz, horrenbeste arduratzeagatik. Gure bizitzetan azaltzeagatik. Peio, beti animoak emateko prest egoteagatik. Nire lorpen txikiekin ni bezain beste pozteagatik. Azken urteetan lanetik kanpoko bizitza errazagoa egiteagatik. Jokin, inoiz aitortu ez arren, txikitatik gauza askotan nire eredu izateagatik. Horrenbeste gauza irakasteagatik. Mundu akademikoan sartzera animatzeagatik. Aizpea, zer esan. Ezagutzen dudan pertsonarik onena izateagatik. Bizitzan ausardia izan behar dela irakasteagatik. Beti nire ondoan egoteagatik. Bihotz-bihotzez, eskerrik asko denoi. Irun-Donostia, 2022ko iraila. Antton Babaze vii LABURPENA Argiaren eta materiaren arteko elkarrekintza aspalditik izan da aztergai zientzialarientzat. Esaterako, 1850eko hamarkadan, Michael Faraday-k aurkitu zuen tamaina nanometrikoko urrezko partikulek eta bolumen handiko urrezko egiturek propietate optiko oso desberdinak dituztela [ 1 ]. Propietate optiko berezi horiek dira, adibidez, urre-errubi beiraren kolore gorri distiratsua sortzen dutenak [ 2 , 3 ]. Aurkikuntza haren azalpen fisikoa Gustav Mie-k eman zuen zenbait urte geroago [ 4 ], 1908an, James Clerk Maxwell-en teoria elektromagnetikoa [ 5 ] erabiliz urrezko nanopartikula txikien sakabanatze-propietateak aztertu zituenean [ 6 – 8 ]. Nanopartikula metaliko txiki hauek argia maiztasun jakin batzuetan (normalean espektro ikusgaian) modu eraginkorrean barreiatzen dutela aurkitu zuen Mie-k; maiztasun hauek materialaren, partikularen tamainaren eta ingurune dielektrikoaren propietateen araberakoak izanik [ 9 ]. Maiztasun jakin horiek gainazaleko plasmoi lokalizatuei dagozkie [ 10 – 12 ], argiaren bidez kitzika daitezkeen nanopartikula metalikoen gainazaleko karga-oszilazioen erresonantziei, alegia [ 13 – 16]. Azken urteetan, nanopartikula metalikoen plasmoi erresonantziek interes handia piztu dute Nanofotonika alorrean [ 17 – 19 ], uhin elektromagnetiko erasotzailearen anplitudea areagotzeko eta argia uhin-luzera baino eskualde txikiagoetan lokalizatzeko duten ahalmena dela-eta [ 20 – 22 ]. Adibidez, bi nanopartikula metalikoren arteko eskualde nanometrikoan (nanobarrunbe plasmoniko deiturikoa, nanogap edota nanocavity ingelesez), eremu elektromagnetikoaren anplitudea 1001000 aldiz handitu daiteke [ 23 ]. Ondorioz, gaur egun plasmoi erresonantziak asko erabiltzen dira hainbat espektroskopia eta mikroskopia tekniketan, hala nola gainazalak areagotutako Raman espektroskopian (surface-enhanced Raman spectroscopy) [ 24 ], gainazalak areagotutako fluoreszentzian (surface-enhanced fluorescence) [ 25 – 27 ], edota molekula bakarren detekzioan (single-molecule imaging) [ 28 , 29 ]. Gainera, plasmoi erresonantziek aplikazio itxaropentsuak dituzte, besteak beste, biomedikuntzan [ 30 – 32 ], energiaren biltegiratzean [ 33 – 35 ], edota optika ez-linealean [36,37]. Efektu plasmonikoak elektromagnetismo klasikoaren teoriaren testuinguruan aztertu izan dira batik bat [ 38 – 45 ], non argiaren eta materiaren arteko elkarrekintza Maxwell-en ekuazioek deskribatzen duten eta sistemaren erantzun optikoa, oro har, linealtzat jotzen den [ 6 ]. Hala ere, egungo konfigurazio esperimentaletan (nanopartikula benetan txikiak, partikulen arteko distantziak azpi-nanometrikoak ix Contents Acknowledgments v Laburpena ix List of abbreviations xv Introduction 1 1 Classical description of light–matter interaction 5 1.1 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.1 The local-response approximation . . . . . . . . . . . . . . . 7 1.1.2 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . 9 1.1.3 The nonretarded approximation . . . . . . . . . . . . . . . . 11 1.2 Plasmonics ............................... 16 1.2.1 Bulkplasmons ......................... 16 1.2.2 Surface plasmons and surface plasmon polaritons . . . . . . 17 1.2.3 Localized surface plasmon polaritons (LSPPs) . . . . . . . . 18 1.3 Nonlinear optical response of small nanostructures . . . . . . . . . 24 1.4 Plexcitonics: Quantum emitter excitons coupled to plasmons . . . 27 1.4.1 The point-dipole approximation . . . . . . . . . . . . . . . . 28 1.4.2 The self-interaction dyadic Green’s function . . . . . . . . . 28 1.4.3 Coupled harmonic-oscillator model . . . . . . . . . . . . . . 31 1.5 Summary ................................ 34 2 Quantum many-body description of light–matter interaction 35 2.1 Fundamentals of density functional theory (DFT) . . . . . . . . . . 36 2.1.1 The local-density approximation (LDA) . . . . . . . . . . . 38 2.1.2 The jellium model of free-electron metals . . . . . . . . . . 39 2.2 Fundamentals of time-dependent density functional theory (TDDFT) 41 2.2.1 The wave-packet propagation (WPP) method . . . . . . . . 42 2.3 Linear optical response of canonical plasmonic nanostructures addressed within TDDFT . . . . . . . . . . . . . . . . . . . . . . . 48 2.3.1 Individual spherical metallic nanoparticles . . . . . . . . . . 49 2.3.2 Dimers of spherical metallic nanoparticles . . . . . . . . . . 55 2.3.3 Cylindrical metallic nanowires . . . . . . . . . . . . . . . . 61 xvi 2.4 Nonlinear effects in the optical response of spherical plasmonic nanoparticles addressed within TDDFT . . . . . . . . . . . . . . . 65 2.5 Semiclassical surface-response formalism (SRF) . . . . . . . . . . . 69 2.5.1 Optical response of spherical nanostructures using the SRF 72 2.6 Summary ................................ 75 3 Quantum surface effects in the electromagnetic coupling between quantum emitters and metallic nanoparticles 77 3.1 Systemandmethods.......................... 78 3.1.1 Time-dependent density functional theory (TDDFT) . . . . 80 3.1.2 Classical local-response approximation (LRA) . . . . . . . . 82 3.1.3 Semiclassical surface-response formalism (SRF) . . . . . . . 82 3.2 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.2.1 Quantum TDDFT vs. classical LRA . . . . . . . . . . . . . 84 3.2.2 Quantum TDDFT vs. semiclassical nondispersive SRF . . . 87 3.2.3 Interpretation of the quantum effects within the nondispersive SRF and its limitations . . . . . . . . . . . . . . . . . . . . 88 3.2.4 Quantum TDDFT vs. nondispersive SRF in a nanoparticle dimer .............................. 91 3.3 Summary ................................ 94 4 Dispersive surface-response formalism to address optical nonlocality in situations of extreme plasmonic field confinement 97 4.1 Systemandmethods.......................... 98 4.2 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.2.1 TDDFT study of the energy dispersion of multipolar plasmon resonances in a metallic nanowire . . . . . . . . . . . . . . . 101 4.2.2 Calculation of the dispersive Feibelman parameter d⊥(ω, k∥)103 4.2.3 Validation of the dispersive SRF . . . . . . . . . . . . . . . 107 4.3 Summary ................................ 112 5 Electronic exciton–plasmon coupling in a nanocavity beyond the electromagnetic interaction picture 115 5.1 Systemandmethods.......................... 116 5.1.1 Characterization of the model quantum emitter (QE) . . . 118 5.2 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.2.1 Influence of the QE exciton quenching at optical frequencies 120 5.2.2 Quantum finite-size effects vs. electronic QE–MNPs coupling 127 5.2.3 Charge-transfer resonances at low frequencies . . . . . . . . 128 5.3 Summary ................................ 130 6 Second-harmonic generation from a quantum emitter coupled to a metallic nanoparticle 131 6.1 Systemandmethods.......................... 132 6.2 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 134 xvii 6.2.1 TDDFTresults......................... 134 6.2.2 TDDFT vs. semi-analytical model . . . . . . . . . . . . . . 138 6.2.3 Influence of the intrinsic losses and the position of the QE on the efficiency of SHG . . . . . . . . . . . . . . . . . . . . 141 6.2.4 Polarization conversion of the second-harmonic field . . . . 146 6.3 Summary ................................ 147 Conclusions and Outlook 149 Appendices A Atomic units 155 B Electric near field induced by a spherical nanoparticle within the classical local-response approximation (LRA) 157 C Dynamics of quantum emitters in the time domain 159 D Dispersive Feibelman parameter obtained for cylindrical nanowires of different size 163 E Plasmon resonances sustained by a cylindrical nanowire within the surface-response formalism 165 F Relationship between the results obtained under finite Gaussian pulses and plane-wave excitation 169 List of publications 173 Bibliography 175 xviii INTRODUCTION The interaction between light and matter has been at the focus of the scientific community over the last centuries. For example, in the 1850s, Michael Faraday discovered that nanoscopic gold colloids exhibit optical properties different to those in bulk metals [ 1 ], giving rise, e.g., to the bright red color of gold-ruby glass [ 2 , 3 ]. Some years later, in 1908, Gustav Mie provided the physical explanation to this effect [ 4 ] when he applied the electromagnetic theory developed by James Clerk Maxwell [ 5 ] to study the scattering properties of small gold nanoparticles [ 6 – 8 ]. Mie found that such small metallic nanoparticles efficiently scatter light at specific frequencies (typically in the visible range) which depend upon the material and size of the nanoparticle as well as upon the properties of the dielectric environment [ 9 ]. It would be later known that these specific frequencies correspond to the so-called localized surface plasmons [ 10 – 12 ], the resonances of the surface charge oscillation sustained by the collective oscillation of the free electrons in metallic nanoparticles that can be excited by light [13–16]. In the last years, the excitation of plasmon resonances in metallic nanoparticles has attracted great interest in Nanophotonics [ 17 – 19 ] due to their capability to enhance and squeeze incident electromagnetic fields into subwavelength regions [ 20 – 22 ]. For example, it is possible to obtain electromagnetic fields up to 1001000 times larger than the incident field in the nanometric region between two metallic nanoparticles (referred to as a plasmonic gap or plasmonic nanocavity) [ 23 ]. As a result, plasmon resonances are widely used in a variety of spectroscopy and microscopy techniques such as surface-enhanced Raman spectroscopy [ 24 ], surface-enhanced fluorescence [ 25 – 27 ], or single-molecule imaging [ 28 , 29 ], and enable promising applications in biomedicine [30–32], energy storage [33–35], and nonlinear optics [36,37], among others. Theoretically, plasmonic effects have been usually studied within the context of classical electromagnetism [ 38 – 45 ], where the interaction between light and matter is described by Maxwell’s equations and the optical response of the system is generally considered to be linear [ 6 ]. However, current experimental configurations involving nanometer-sized metallic nanoparticles and interparticle distances push light–matter interaction to the limit where quantum many-body phenomena influence optical properties [ 46 – 52 ]. In these extreme situations, classical descriptions are no longer valid [ 53 – 55 ], and alternative theoretical approaches allowing for incorporating nonlocality [ 56 – 60 ], electronic spill in/out [ 61 – 63 ], surface-enabled Landau damping 1 Introduction [ 64 – 66 ], and electron tunneling [ 67 , 68 ] are required to describe the optical properties of miniaturized plasmonic structures correctly. In this context, time-dependent density functional theory (TDDFT) [ 69 – 72 ] is the method of choice to address the optical response of nanoscale systems in this thesis, since this theory accounts for the quantum nature of the electron dynamics from first principles by addressing the time evolution of the electron density in metallic nanoparticles subjected to external illumination [ 73 ]. In addition, the use of TDDFT calculations is not restricted to obtaining only the linear optical response, but it also enables to directly obtain the nonlinear response of plasmonic systems [ 74 – 78 ]. TDDFT is thus a very powerful tool to capture quantum many-body phenomena involved in the optical and electronic (i.e., optoelectronic) response, and it sets up the cornerstone of this thesis. On the other hand, since TDDFT is limited to addressing small systems that contain a few thousands of atoms due to its computational complexity [79,80], less-demanding semiclassical models have also been developed to capture various quantum effects [ 81 , 82 ]. Here we can mention, for instance, the quantum-corrected model accounting for electron tunneling in subnanometric metallic junctions [ 83 – 87 ], several hydrodynamic descriptions incorporating nonlocality of the optical response of metallic nanostructures [ 88 – 94 ], and the surface-response formalism (SRF) [ 95 – 97 ] based on the inclusion of quantum surface-response corrections (the so-called Feibelman parameters) at the metal–dielectric boundaries [ 98 ]. Together with TDDFT, in this thesis we also perform simulations based on the SRF, which captures important quantum effects such as the spill in/out of the induced charge density and surface-enabled Landau damping, but cannot account e.g. for electron tunneling between two nanoparticles in close proximity. TDDFT and the aforementioned semiclassical models are often used to study the optical response of isolated plasmonic nanostructures, but they can also be used to analyze quantum effects in the optoelectronic interaction between plasmonic nanostructures and quantum emitters such as atoms, quantum dots or organic molecules. Indeed, the optical response of a coupled emitter–plasmonic nanostructure system has been widely studied in Nanophotonics. The plasmonic field enhancement and confinement allow quantum emitters to interact much more efficiently with light, leading to a plethora of interesting phenomena [ 99 , 100 ]. Some of these effects can be described by classical calculations, such as the enhancement of the spontaneous emission rate (Purcell effect [ 101 , 102 ]) and the modification of the resonant frequency (Lamb shift [ 103 , 104 ]) of a quantum emitter located near a plasmonic nanostructure. On the other hand, the coupling between a plasmonic cavity and a quantum emitter can also involve quantummechanical interactions that can affect the chemical properties of the emitter [ 105 ], induce electron transfer between the nanoparticles and the emitter [ 106 – 109 ], and eventually modify drastically the nonlinear optical response of the coupled system [110,111]. This thesis aims at theoretically studying novel quantum many-body phenomena in the optoelectronic response of plasmonic nanostructures and their interaction with quantum emitters. In particular, we seek a deeper fundamental knowledge 2 Introduction into quantum-mechanical effects that occur in plasmonic cavities of nanometric size or smaller, where nonlocality, nonlinearity, or electron-transfer effects can all play an important role. To this end, quantum TDDFT simulations are used primarily, but classical calculations based on the local-response approximation within the dielectric framework as well as semiclassical models based on the SRF are also employed to complete the analysis and provide additional insights. In Chapter 1, we review the fundamentals of classical light–matter interaction and explain the key concepts of this thesis from a classical perspective, focusing on the classical description of localized surface plasmon polaritons sustained by spherical metallic nanoparticles and their dimers as well as on the electromagnetic interaction between plasmonic nanostructures and quantum emitters. In Chapter 2, we describe the quantum TDDFT methodology and the semiclassical SRF adopted in this thesis to account for quantum many-body effects present in the light–matter interaction, and introduce quantum-mechanical concepts such as electron tunneling, surface-enabled Landau damping, and electron spill out with the examples of the optical response of canonical plasmonic nanostructures. In Chapter 3, we use TDDFT to study the influence of quantum surface effects in the electromagnetic interaction between quantum emitters and metallic nanoparticles. We identify the dominant effects when electron tunneling and the electronic coupling between the quantum emitter and the metallic nanoparticle are negligible. We further analyze the origin of the observed effects with semiclassical calculations based on a SRF that neglects the optical nonlocality in the direction parallel to the metallic surface, as usually implemented in the literature, and establish the validity range of this standard implementation of the semiclassical SRF. In Chapter 4, we show that the standard implementation of the SRF can be extended by explicitly accounting for the nonlocality of the optical response in the parallel direction along the metallic surface in the calculation of the Feibelman parameters. The inclusion of this nonlocality has not been so far considered in the literature to the best of our knowledge, and produces a broader range of validity of the SRF including extreme subnanometric configurations. In Chapter 5, we analyze the effect of electronic coupling between the electronic states of a quantum emitter and those of a metallic nanoparticle dimer when the distance between the nanoparticles and the emitter is subnanometric. This effect has been neglected in Chapters 3and 4, which focused on the electromagnetic interaction. We demonstrate that, in such situations, the electronic interaction between the emitter and the nanoparticles drastically modifies the optical resonances of the coupled system. Finally, in Chapter 6, we study the nonlinear optical response of a system consisting of a quantum emitter and a spherical metallic nanoparticle, and show that the electromagnetic emitter– nanoparticle coupling can enable strong nonlinear second-harmonic generation, otherwise forbidden due to symmetry constrains if an isolated spherical nanoparticle were considered. The content of this thesis thus presents a quantum mechanical many-body approach to the optical response of plasmonic cavities, which provides new insights into coupling with emitters, electron transfer processes and nonlinear effects. All these effects are of paramount importance in nowadays state-of-the-art 3 Introduction Nanophotonics. 4 Chapter 1 CLASSICAL DESCRIPTION OF LIGHT–MATTER INTERACTION 1.1 Maxwell’s equations The interaction between light and matter is well addressed from a classical and macroscopic point of view by the time-dependent Maxwell’s equations [6,22], ∇×E(r, t) = −∂ ∂tB(r, t),(1.1a) ∇×H(r, t) = ∂ ∂tD(r, t) + Jext(r, t),(1.1b) ∇·D(r, t)=4πρext(r, t),(1.1c) ∇·B(r, t)=0,(1.1d) which represent a set of coupled partial differential equations where light is described as an electromagnetic wave with electric E ( r, t )and magnetic H ( r, t )field components. 1 The time-dependent Maxwell’s equations in Eq. (1.1) determine the dynamics of the electromagnetic wave (light) at a position r and instant of time t through a given dielectric material in response to an external charge density ρext ( r, t ) and current density Jext ( r, t )that act as sources of electromagnetic radiation. For the nonmagnetic materials that we consider in this thesis, the magnetic induction B ( r, t )is straightforwardly linked to H ( r, t )by B ( r, t ) = 1 c2H ( r, t ), with c the speed of light in a vacuum. Finally, the displacement vector, D ( r, t ), is related to the electric field E ( r, t )within the linear-response regime by the following constitutive 1 Atomic units (au) are used throughout this manuscript unless otherwise stated (see Appendix A). 5 Chapter 1. Classical description of light–matter interaction relationship D(r, t) = Z∞ −∞ Z∞ −∞ dr′dt′ε(r−r′, t −t′)E(r′, t′),(1.2) where ε ( r−r′, t −t′ )is the dielectric function describing optical excitations in the material. We assume in Eq. (1.2) that the medium is isotropic and homogeneous by considering that the dielectric function ε is a scalar quantity that spatially depends on r−r′ . Moreover, the causality of the dielectric response is introduced by imposing ε(r−r′, t −t′)=0for any t′> t. According to Eq. (1.2) , the displacement vector D at a particular position r and instant of time t depends on the value of the electric field E at all positions r′ and times t′< t . Thus, the dielectric response of a material is, in general, nonlocal both in space and time. To deal with the temporal nonlocality, it is convenient to use the time-to-frequency Fourier transform 2 and express Eq. (1.2) in the domain of the angular frequency ωof the electromagnetic field, D(r, ω) = Z+∞ −∞ dr′ε(r−r′, ω)E(r′, ω).(1.4) Consequently, by applying the same Fourier transform to Eq. (1.1) , Maxwell’s equations can be expressed in the frequency domain as: ∇×E(r, ω) = iωB(r, ω),(1.5a) ∇×H(r, ω) = −iωD(r, ω) + Jext(r, ω),(1.5b) ∇·D(r, ω)=4πρext(r, ω),(1.5c) ∇·B(r, ω)=0.(1.5d) Equations (1.1) and (1.5) are, of course, equivalent, and can be solved by adopting different techniques. Among numerical methods to solve Maxwell’s equations, we can cite the Boundary Element Method [ 112 ], the Finite-Difference Time-Domain method [ 113 ], the Discrete-Dipole Approximation [ 114 ] or the Finite-Element Method [115]. On the other hand, the influence of the spatial nonlocality of the dielectric function present in Eq. (1.4) on the optical response of different materials is currently an active research topic by itself [ 81 ], and many efforts have been devoted 2 In this thesis, we use the following definitions for the time-to-frequency ( F ) and for the frequency-to-time (F−1) Fourier transforms of a function f: f(ω) = F[f(t)] = Z+∞ −∞ f(t)eiωt dt, f(t) = F−1[f(ω)] = 1 2πZ+∞ −∞ f(ω)e−iωt dω, (1.3) 6 1.1. Maxwell’s equations to the development of numerical methods that correctly account for such nonlocality [ 59 , 60 , 88 , 89 , 91 , 92 , 94 , 116 – 123 ]. Using these methods, it has been shown that spatial nonlocality can play a significant role in defining the properties of the optical response of metallic nanostructures with characteristic dimensions below ∼ 10 nm [ 56 , 86 , 124 ], as it is the case of the systems studied in this thesis. Indeed, we discuss in Chapters 2,3,4and 5the role that spatial nonlocality (and other quantum phenomena) plays in different scenarios involving small metallic nanoparticles when they interact with light and with quantum emitters. We account for such effects by using quantum time-dependent density functional theory (TDDFT) simulations (see Chapter 2), and, at a different level of approximation, by a semiclassical surfaceresponse formalism (SRF) incorporating the Feibelman parameters obtained ab initio (Section 2.5). In contrast, the classical calculations performed within this thesis are based on a local description of the optical response of metals, so that nonlocal effects can only be partially introduced with the use of phenomenological parameters [58,125,126], as we detail in the following Subsection 1.1.1. 1.1.1 The local-response approximation The local-response approximation (LRA) neglects the spatial nonlocality of the dielectric function ε ( r−r′, ω )by assuming that the dielectric response of a material at a position r is independent of the response at any other position r′ = r . This assumption can be expressed as ε(r−r′, ω) = ε(r, ω)δ(r−r′),(1.6) with δ ( r−r′ )the Dirac delta function. The displacement vector D ( r, ω )thus transforms from Eq. (1.4) to a much simpler constitutive relationship, D(r, ω) = Z+∞ −∞ dr′ε(r, ω)δ(r−r′)E(r′, ω) =ε(r, ω)E(r, ω). (1.7) Despite its simplicity, the LRA has successfully described many physical phenomena of interest in the context of light–matter interaction at the nanoscale, particularly in situations where the characteristic length scale of the studied nanostructures is substantially larger than the Fermi wavelength of electrons in metals. The Drude response model of free-electron metals The simplest way to estimate the (frequency-dependent) local dielectric function ε ( r, ω )of a metal is to adopt the Drude model [ 127 , 128 ]. This model considers that conduction electrons freely move within a homogeneous gas in response to an external electric field E ( t )while the heavy metal ions remain immobile. The equation of motion of a conduction electron is then ¨ x(t) + γp˙ x(t) = −E(t),(1.8) 7 Chapter 1. Classical description of light–matter interaction Figure 1.2:(a) Spherical coordinate system ( r, θ, φ )with the origin at the center of a spherical nanoparticle of radius a characterized by a dielectric function ε ( ω ).(b) Coordinate system employed to address the optical response of a spherical dimer formed by two nanoparticles of radius a1 and a2 separated by a gap of size D . The nanoparticles are characterized with a dielectric function ε1 ( ω )and ε1 ( ω ). The coordinates are expressed in spherical coordinates ( r1, θ1, φ1 )and ( r2, θ2, φ2 )with the origin at the center of each nanoparticle. In both individual nanoparticle and dimer structures, the system is surrounded by a dielectric material characterized with ϵd. Since a point dipole of amplitude p ( ω )located at the center of the coordinate system creates an electrostatic potential at a position r given by ϕind ( r, ω ) = r·p(ω) |r|3 , p(ω)can be obtained as p(ω) = m=1 X m=−1 b1m(ω)a2Ym 1(π 2,0) ˆ x+Ym 1(π 2,π 2)ˆ y+Ym 1(0,0) ˆ z,(1.29) where {ˆ x,ˆ y,ˆ z}are the unit vectors along the {x, y, z}-axes. Optical response of a spherical dimer We now consider the case of a dimer consisting of two spherical nanoparticles of radius a1 and a2 separated by a gap of size D (see Figure 1.2b). Each nanoparticle is characterized by a dielectric function ε1 ( ω )and ε2 ( ω ), and the entire system is surrounded by a dielectric material characterized by εd . The numerical implementation described in this section to obtain the optical response of this system is based on a coupled-multipole method [ 131 , 132 ], which uses the solution of the individual spherical nanoparticle explained above. The potential ϕind dimer ( r, ω )induced by the dimer in response to an external potential ϕext ( r, ω )is given by the sum of the potential induced by each nanoparticle 14 1.1. Maxwell’s equations (1 and 2), ϕind dimer(r, ω) = ϕind 1(r1, ω) + ϕind 2(r2, ω),(1.30) where the vectors r1≡ ( r1sinθ1cosφ1, r1sinθ1sin φ1, r1cosθ1 )and r2≡ ( r2sinθ2cosφ2, r2sinθ2sin φ2, r2cosθ2 )are written in spherical coordinates with the origins at the center of the corresponding nanoparticle, as shown in Figure 1.2b. We define the coordinates such that the centers of the nanoparticles are located at the z -axis separated by a distance δ = a1 + a2 + D . Thus, r1 and r2 are related by r2=qr2 1+δ2−2δr1cosθ1, cosθ2= (r1cosθ1−δ)/r2, φ2=φ1. (1.31) The electrostatic potential ϕind i ( ri, ω )induced by the nanoparticle i (with i = 1 , 2) follows identical expression as that of the individual nanoparticle explained above [Eq. (1.26)], ϕind 1(r1, ω) = ℓmax X ℓ=0 ℓ X m=−ℓ bℓm 1(ω)Ym ℓ(θ1, φ1)rℓ 1/aℓ 1r1< a1 aℓ+1 1/rℓ+1 1r1> a1,, ϕind 2(r2, ω) = ℓmax X ℓ=0 ℓ X m=−ℓ bℓm 2(ω)Ym ℓ(θ2, φ2)rℓ 2/aℓ 2r2< a2 aℓ+1 2/rℓ+1 2r2> a2,, (1.32) with, in this case, bℓm 1(ω) = ξ1(ω, ℓ)Zd˜ Ω1[Ym ℓ(θ1, φ1)]∗∂ ∂r1ϕext(r1, ω) + ϕind 2(r2, ω)r1=a1 , bℓm 2(ω) = ξ2(ω, ℓ)Zd˜ Ω2[Ym ℓ(θ2, φ2)]∗∂ ∂r2ϕext(r2, ω) + ϕind 1(r1, ω)r2=a2 , (1.33) where ξ1(ω, ℓ) = −a1(ε1−εd) εd(ℓ+ 1) + ε1ℓ, ξ2(ω, ℓ) = −a2(ε2−εd) εd(ℓ+ 1) + ε2ℓ. (1.34) The integrals in Eq. (1.33) take into account that the potential induced by one nanoparticle acts as external potential for the other one. From Eq. (1.32) , the coefficients bℓm 1 ( ω )and bℓm 2 ( ω )in Eq. (1.33) can be found from the following expression given in matrix form: bℓm 1(ω) = I−T2→1T1→2−1bℓm,ext 1(ω) + T2→1bℓm,ext 2(ω), bℓm 2(ω) = I−T1→2T2→1−1bℓm,ext 2(ω) + T1→2bℓm,ext 1(ω). (1.35) 15 Chapter 1. Classical description of light–matter interaction Here, bℓm 1 ( ω )and bℓm 2 ( ω )represent column vectors containing the coefficients bℓm 1 ( ω ) and bℓm 2 ( ω ), I is the identity matrix, and bℓm,ext 1 and bℓm,ext 2 are column vectors with the following elements: bℓm,ext 1(ω) = ξ1(ω, ℓ)Zd˜ Ω1[Ym ℓ(θ1, φ1)]∗∂ ∂r1 ϕext(r1, ω)r1=a1 , bℓm,ext 2(ω) = ξ2(ω, ℓ)Zd˜ Ω2[Ym ℓ(θ2, φ2)]∗∂ ∂r2 ϕext(r2, ω)r2=a2 . (1.36) Finally, the elements (ℓm, ℓ′m′)of matrices T2→1and T1→2are given by T2→1 ℓm,ℓ′m′=ξ1(ω, ℓ)Zd˜ Ω1[Ym ℓ(θ1, φ1)]∗∂ ∂r1 Ym′ ℓ′(θ2, φ2)aℓ′+1 2 rℓ′+1 2!r1=a1 , T1→2 ℓm,ℓ′m′=ξ2(ω, ℓ)Zd˜ Ω2[Ym ℓ(θ2, φ2)]∗∂ ∂r2 Ym′ ℓ′(θ1, φ1)aℓ′+1 1 rℓ′+1 1!r2=a2 . (1.37) Once the coefficients bℓm 1 ( ω )and bℓm 2 ( ω )are obtained by solving Eq. (1.35) , the physical quantities of interest such as the induced electric field Eind ( r, ω )or the dipole moment p ( ω )of the dimer structure can be determined as the sum of the corresponding contributions of each nanoparticle [Eqs. (1.28) and (1.29)]. 1.2 Plasmonics Plasmons are collective oscillations of the free-electron gas inside a metal [ 22 , 133 ]. From a classical electromagnetic point of view, plasmons manifest themselves as absorption and scattering resonances obtained from the solutions of Maxwell’s equations and the associated boundary conditions discussed in Section 1.1. On the other hand, within quantum-electrodynamics, the word plasmon refers to the quantum of the metallic electron plasma oscillation [134,135]. Depending on the geometry of the system, plasmons can be excited in the bulk or at the metal–dielectric boundary, and, in the latter situation, they can be either propagating or localized electron-density oscillations. In the following, we describe some canonical examples of plasmonic excitations within the classical electromagnetic theory. In addition to providing the fundamental concepts of plasmonics, this classical description establishes a reference point to compare with the results obtained within a quantum many-body approach, thus allowing us to identify the origin of the effects of interest in this thesis. 1.2.1 Bulk plasmons Bulk or volume plasmons are longitudinal excitations consisting in the coherent oscillation of the electron gas propagating in an infinitely extended metal, where the excited electrons move collectively oscillating at the same frequency. These 16 1.2. Plasmonics longitudinal electromagnetic waves ( ∇·E ( r, ω )  = 0) can be described within the local-response approximation (see Subsection 1.1.1), and can only exist if ε ( ω ) = 0 [ 54 , 57 ]. Thus, assuming a lossless Drude dielectric function [Eq. (1.10) , with γp = 0], the bulk plasmon frequency is equal to the plasma frequency, ωbulk = ωp [Eq. (1.12)]. 1.2.2 Surface plasmons and surface plasmon polaritons When we consider a semi-infinite planar metal slab ( ε ( r, ω ) = ε ( ω ) , z < 0) in contact with vacuum ( ε ( r, ω ) = εd = 1 , z > 0), the translational invariance of the bulk material along the z -direction is broken. As a consequence, a new type of plasmon emerges associated with an oscillating charge density propagating along the metal–dielectric interface. These density oscillations are referred to as surface plasmons. The resonance condition of a surface plasmon can be derived within the nonretarded approximation from the self-sustained induced potential ϕind ( r, ω ) [Eq. (1.22), with ϕext(r, ω)=0]. ϕind(r, ω)is given in this case by [10] ϕind(r, ω) =    ϕ<(ω)eik∥·r∥e|k∥|zz < 0 ϕ>(ω)eik∥·r∥e−|k∥|zz > 0 ,(1.38) where r∥ and k∥ are the two-dimensional position and momentum vectors in the plane of the metal surface, respectively, and ϕ< ( ω )and ϕ> ( ω )are frequencydependent coefficients to be determined from the boundary conditions given by Eq. (1.23) . From the continuity of the potential at z = 0 [Eq. (1.23a) ] it is obtained that ϕ< ( ω ) = ϕ> ( ω ). Further, from the continuity of the normal component of D ( r, ω ) = −ε ( r, ω ) ∇ϕ ( r, ω )at z = 0 [Eq. (1.23b) ], the resonance condition of the nonretarded surface plasmon is given by ε ( ω ) = − 1. Thus, assuming a lossless Drude dielectric function [Eq. (1.10) , with γp = 0], the long-wavelength limit of the surface plasmon frequency is equal to ωsp =ωp/√2.(1.39) On the other hand, a surface plasmon polariton (SPP) is a hybrid mode that propagates at the metal–dielectric boundary and that results from the coupling of a surface plasmon and an electromagnetic wave [ 136 ]. These SPPs are characterized by electromagnetic fields that are evanescent in the direction normal to the surface, thus exponentially decaying with increasing distance from the interface (see sketch in Figure 1.3a). Considering an electromagnetic wave that satisfies Maxwell’s equations in this system [Eq. (1.5) ], and applying the boundary conditions given by Eqs. (1.16) and (1.18) as well as the conservation of momentum vector k∥ along the propagation direction [ 22 ], the following dispersion relationship of SPPs is obtained (for a lossless Drude dielectric function), ω2 spp(k∥) = ω2 p/2 + c2k2 ∥−qω4 p/4 + c4k4 ∥,(1.40) 17 Chapter 1. Classical description of light–matter interaction + + + − − − Figure 1.3:(a) Schematic representation of a SPP with a wavenumber kspp propagating along the metal–vacuum interface located at z = 0. The electromagnetic field E associated with the SPP is represented by brown lines. (b) Frequency of a SPP obtained in this system from Eq. (1.40) ( ωspp , blue line), of light in free space ( ωc , red line), and of a long-wavelength nonretarded surface plasmon ( ωsp , green line) as a function of the wavenumber k∥ parallel to the metal–vacuum interface. The frequency ω in the vertical axis is measured in units of ωp , and k∥ in the horizontal axis in units of ωp/c, with cthe speed of light in vacuum. with c the speed of light in a vacuum. The amplitude of the momentum vector k∥ in Eq. (1.40) corresponds to the wavenumber kspp of a SPP, kspp = k∥ , as schematically depicted in Figure 1.3a. As shown in Figure 1.3b, for small k∥ , the SPP frequency ωspp (blue line) approaches the dispersion line of light in free space, ωc = ck∥ (red line). In contrast, for large k∥ , ωspp yields the nondispersive ( k∥ -independent) frequency of the nonretarded surface-plasmon frequency in the long-wavelength limit [Eq. (1.39) ], ωsp = ωp/√2 (green line) [ 10 ]. Importantly, for a given frequency ω , the momentum k∥ of light in free space is always smaller than the momentum kspp of SPPs, and therefore, a laser beam incident on an ideal surface cannot excite SPPs because momentum and energy cannot be simultaneously preserved. For this reason, several mechanisms have been adopted to provide the extra momentum needed to excite SPPs with light, such as the use of surface roughness or gratings [ 137 ], evanescent fields [ 138 , 139 ], or sharp metallic tips placed on top of the metal–dielectric surfaces [140]. 1.2.3 Localized surface plasmon polaritons (LSPPs) When metallic nanoparticles (MNPs) of finite size are considered, a new type of plasmon resonances emerges, so-called localized surface plasmon polaritons (LSPPs) [ 11 , 141 ]. The word localized is used because, unlike the SPP previously discussed, in this case the surface-charge oscillations are not propagating in space. Indeed, the transnational invariance of the system is broken, and plasmons are confined in the three dimensions of space. Thus, LSPPs can often be understood as confined SPPs with a quantized wavenumber kspp . Importantly, contrary to SPPs in semi-infinite 18 1.2. Plasmonics metal slabs, LSPPs in small MNPs can be excited by an incident laser beam, because the finite geometry of the MNP provides the extra momentum needed to overcome the mismatch between the momentum of light and plasmons. For simplicity, in this thesis we refer to LSPPs as “localized surface plasmons (LSPs)”, or simply “plasmon resonances” or “plasmons”. In this thesis, we focus on quantum effects associated with the excitation of LSPPs in small MNPs, where special attention is paid to spherical nanoparticles and nanoparticle dimers. In the following, we describe the basic properties of LSPPs sustained in such geometries using a classical electromagnetic theory. LSPPs in spherical nanoparticles We first derive the resonance conditions fulfilled by LSPPs in a spherical MNP of radius a surrounded by vacuum (dielectric function εd = 1). These conditions are again determined from the self-sustained induced potential ϕind ( r, ω )that satisfies Laplace’s equation [Eq. (1.22) ], which for the spherical MNP is given by Eq. (1.26) . Under no external excitation, ϕext ( r, ω )=0, the coefficients bℓm ( ω )are non-null only at the poles of Eq. (1.27) , so that the resonance conditions of the nonretarded LSPPs sustained in a spherical MNP are given by ℓε(ω)+(ℓ+ 1) = 0,(1.41) where ℓ is the multipole order of the LSPP resonance. Therefore, assuming a lossless Drude dielectric function [Eq. (1.10) , with γp = 0], the frequencies ωℓ of LSPPs in a spherical MNP are given by: ωℓ=ωprℓ 2ℓ+ 1,(1.42) which are independent of the size of the nanoparticle. However, this expression is rigorous only within the validity range of the nonretarded and local-response approximations. Indeed, it is known that the scale invariance of Eq. (1.42) is lifted when retardation effects are considered in large MNPs ( a≳ 15 nm ) [ 7 , 142 ], as well as when quantum-size effects are considered in very small nanostructures ( a≲ 5 nm ) [ 61 , 95 , 96 ]. We return in Chapters 2,3, and 4to the impact of quantum-size effects on the resonant frequencies of LSPPs, ωℓ. Figure 1.4a shows the values of the LSPP frequencies ωℓ as a function of the multipole order ℓ [Eq. (1.42) ]. For small values of ℓ , the resonance frequencies ωℓ of different order ℓ can be well differentiated. For example, the dipolar plasmon (DP) frequency ( ℓ = 1) is given by ωDP = ωℓ=1 = ωp/√3 , the quadrupolar plasmon (QP) frequency ( ℓ = 2) by ωQP = ωℓ=2 = ωpq2 5 , and the octupolar plasmon (OP) frequency ( ℓ = 3) by ωOP = ωℓ=3 = ωpq3 7 . In contrast, different LSPPs with large ℓ have very similar resonant frequencies, and thus, due to the broadening of each resonance given by the losses ( γp in a Drude metal), they pile up in a single broad resonance, the so-called pseudomode [ 143 , 144 ] (see Chapter 3). The 19 Chapter 1. Classical description of light–matter interaction + + + − − − +− Figure 1.4:(a) Frequencies ωℓ of the nonretarded LSPP resonances sustained in a spherical MNP as a function of the multipole order ℓ , considering a lossless Drude dielectric function [Eq. (1.10) ]. The frequencies ωℓ in the vertical axis are given in units of ωp . The green dashed line represents the frequency of the surface plasmon, ωsp/ωp = 1 /√2 .(b) Sketch of the surface charge density induced in a spherical MNP by an external electromagnetic field Eext corresponding to plane-wave illumination. The green arrow represents the oscillating dipole moment p induced at the MNP. (c) Absolute value of the electric field Eind induced by a spherical MNP in response to an external plane wave Eext , calculated within the nonretarded approximation. The MNP has a radius a = 5 nm , and it is characterized by a Drude dielectric function with ωp = 5 . 89 eV and γp = 0 . 21 eV . The color map is shown at the DP resonance, ωDP = 3 . 4 eV .(d) Extinction (blue) and scattering (red, multiplied by 10) cross-section spectra of the same MNP as in panel (c). Due to the small size of the MNP, the extinction and absorption cross sections are almost equivalent, σext(ω)≈σabs(ω). frequency of this pseudomode approaches the frequency of the surface plasmon, ωℓ→∞ ∼ωp/√2 = ωsp [Eq. (1.39) ], since for large ℓ the local curvature of the spherical MNP is very large (approximately flat interface) as compared to the wavelength of the excited LSPPs, λLSPP ∼2πa/ℓ. We next show that a plane wave illuminating a spherical MNP can excite LSPPs. We consider a plane wave linearly polarized along the z -direction, which is expressed within the nonretarded approximation as Eext ( r, ω ) = E0ˆ z . Intuitively, Eext ( r, ω )polarizes the MNP along the z -direction, producing a displacement of the free-electron gas with respect to the positively charged background (see sketch in Figure 1.4b). Then, seeking to restore the equilibrium, the Coulomb interaction between the positive and negative charge densities produces a collective oscillation of the free-electron gas. The LSPP corresponds to the resonant excitation of this 20 1.2. Plasmonics collective oscillation. A quantitative study of the excitation of LSPPs under plane-wave illumination can be done by calculating the optical response of the MNP following the procedure explained in Subsection 1.1.3. The external potential ϕext ( r, ω )in this case is given in spherical coordinates by: ϕext(r, ω) = −E0rcosθ, (1.43) and the integral in Eq. (1.27) providing the coefficients bℓm ( ω )has an analytical solution [ 145 ], which is non-null only for ℓ = 1. This means that, in contrast with SPPs, the DP resonance can be excited by plane-wave illumination. The electric field Eind(r, ω)induced in the proximity of the MNP is given by: Eind(r, ω) = −∇ϕind(r, ω) =          −E0ε(ω)−1 ε(ω)+2 (cosθˆ r−sinθˆ θ) | {z } ˆ z r < a E0ε(ω)−1 ε(ω)+2 a3/r3(2 cosθˆ r+ sin θˆ θ)r > a , (1.44) with ˆ r and ˆ θ the unit vectors along the radial and tangential directions, respectively. Note that, as expected, Eind ( r, ω )possesses a resonance at the DP frequency, ωDP =ωp/√3. We show in Figure 1.4c the absolute value of the electric field Eind ( r, ω )induced at the DP frequency, ωDP = 3 . 4 eV , in the proximity of a MNP of radius a = 5 nm . The MNP is characterized by a Drude dielectric function [Eq. (1.10) ] using the parameters ωp = 5 . 89 eV and γp = 0 . 21 eV to describe sodium. Results are normalized to the amplitude of the incoming field, E0 . Crucially, the induced field is much stronger than the incident field, and it is localized in a space region much smaller ( ∼ 1 − 5 nm ) than the wavelength of the incident field ( ∼ 400 nm ). Thus, this electromagnetic-field localization surpasses the diffraction limit. Moreover, we note that the electric field Eind ( r, ω )induced outside the MNP boundary ( r > a ) corresponds to the electrostatic field created by a point dipole placed at the center of the MNP [ 6 , 22 ]. Indeed, the external field Eext induces a dipole moment p ( ω ) = α ( ω ) Eext at the spherical MNP that depends on the quasi-static polarizability [22] α(ω) = a3ε(ω)−1 ε(ω)+2.(1.45) This induced dipole p ( ω )emits light into the far field. The MNP thus acts as an optical nanoantenna that can extremely localize incident electromagnetic radiation in the near field and radiate it into the far field [ 146 ]. The power Psca scattered into the far field is related to p(ω)by [6] Psca =ω4 3c3|p(ω)|2.(1.46) 21 Chapter 1. Classical description of light–matter interaction Typically, the scattered power Psca = σscaI0 is normalized to the intensity I0 = c|Eext|2/ (8 π )of the incident plane wave, leading to the scattering cross section [7], σsca(ω) = 8π 3ω c4|α(ω)|2.(1.47) Moreover, the total power provided by the incoming light to a MNP, Pext , is given by the extinction cross section σext through the relationship Pext = σextI0 . σext(ω)is defined as [7] σext(ω) = 4πω cIm{α(ω)}=σabs(ω) + σsca(ω),(1.48) and corresponds to the sum of both scattering σsca and absorption σabs cross sections, the latter determining the power of light absorbed by the MNP, Pabs =σabsI0. Figure 1.4d displays the extinction σext ( ω )(blue line) and scattering σsca ( ω ) (red line) cross-section spectra of the same spherical sodium MNP considered in Figure 1.4c ( a = 5 nm ). For such small MNP, σsca ( ω ) ∝a6ω c4 (with ωa c<< 1) is about two-orders of magnitude smaller than σext ( ω ) ∝a3ω c . As a consequence, the extinction of the MNP is completely governed in this case by its absorption. Thus, for the small size of the MNPs considered in this thesis, σabs ( ω )and σext ( ω ) are almost equivalent. We therefore use the following expression to compute the absorption cross section in all the thesis [see Eq. (1.48)]: σabs(ω) = 4πω cIm{α(ω)}.(1.49) LSPPs in nanoparticle dimers When two spherical MNPs are placed in close proximity, they form a dimer sustaining coupled LSPPs that can be excited by plane-wave illumination [ 23 ]. The study of LSPPs in the dimer configuration is particularly interesting because they induce much stronger field enhancements than individual MNPs [ 147 , 148 ]. Moreover, the analysis of the optical properties of a spherical dimer can help to understand the physics behind the optical response of more complex plasmonic nanostructures. The resonant LSPP frequencies in a metallic dimer differ from those of the individual MNPs because the electromagnetic coupling between the two MNPs modifies the resulting plasmonic resonances [ 149 ]. The LSPP frequencies in this dimer configuration depend on the material and size of the MNPs forming the dimer, as well as the surface-to-surface gap separation D . As an example, we show in Figure 1.5a the waterfall spectra of the absorption cross section, σabs ( ω ) [Eq. (1.49) ], of a MNP dimer with gap separation ranging from D = 1 . 05 nm (bottom) to D = 3 . 6 nm (top). We consider the same sodium MNPs as in Figures 1.4c-d, and the external illumination is polarized along the dimer axis. Results are 22 1.2. Plasmonics 0 500 1000 1500 2000 2500 3000 Figure 1.5:(a) Waterfall spectra of the absorption cross section σabs ( ω )of a dimer composed by two spherical sodium MNPs of radius a = 5 nm characterized by a Drude dielectric function with ωp = 5 . 89 eV and γp = 0 . 21 eV (same as in Figures 1.4c-d). The gap separation is varied from D = 1 . 05 nm (bottom) to D = 3 . 6 nm (top) in steps of ∼ 0 . 21 nm .(b) Sketch of the optical hybridization between LSPPs of individual MNPs leading to the bonding dipolar plasmon (BDP) and the bonding quadrupolar plasmon (BQP). (c) Color map of the induced electric-field enhancement |E|ind/E0 in the middle of the gap formed by two MNPs with radius a = 5 nm as a function of the frequency of the external plane-wave illumination, ω , and the surface-to-surface gap distance, D . Results are obtained under the nonretarded approximation considering an external electromagnetic field, Eext =E0ˆ z, polarized along the dimer axis (z-axis). obtained within the nonretarded approximation following the procedure explained in Subsection 1.1.3. Two distinct resonances can be observed in the absorption spectra of Figure 1.5a for the whole range of gap separation D considered: an intense bonding dipolar plasmon (BDP) shifting from ωBDP ∼ 3 . 2 eV to ∼ 2 . 8 eV as D is reduced, and a weaker bonding quadrupolar plasmon (BQP) at ωBQP ∼ 3 . 5 eV , more pronounced for narrow gaps. The BQP also redshifts with decreasing D . The origin of these plasmon resonances can be understood using a hybridization picture that considers the coupling between the multipole modes ℓ of the two individual MNPs [Eq. (1.42) ] similarly to the hybridization between atomic orbitals in diatomic molecules [ 150 ]. According to this hybridization picture (see sketch in Figure 1.5b), the BDP is mainly created due to the electromagnetic coupling between the DP resonances ( ℓ = 1) of the individual MNPs, while the BQP is mainly a consequence of the mixture between the DP mode of one MNP and the quadrupolar mode ( ℓ = 2) of 23 Chapter 1. Classical description of light–matter interaction + − Figure 1.7:(a) Absorption cross-section spectra σabs ( ω )of an individual spherical MNP (dashed black line) and that of a QE–MNP coupled system (blue line). The MNP has a radius a = 5 nm , and it is characterized by a Drude dielectric function with ωp = 5 . 89 eV and γp = 0 . 21 eV . The QE, located at 1 . 5 nm from the surface of the MNP, is characterized by an oscillator strength α0 QE = 2 au , an intrinsic loss rate γQE = 10 meV , and a transition frequency ωQE = ωDP = 3 . 4 eV in resonance with the DP of the MNP. (b) Same as in panel (a) but for a MNP dimer characterized by a gap separation D = 3 nm . The QE emitter is at resonance with the BDP mode of the dimer, ωQE =ωBDP = 3.15 eV. The results are obtained within the nonretarded approximation. by the corresponding components [145], G⊥(rQE,rQE, ω) = ∞ X ℓ=1 ε(ω)−1 ε(ω) + ℓ+1 ℓ a2ℓ+1 R2ℓ+4 (ℓ+ 1)2,(1.58a) G∥(rQE,rQE, ω) = ∞ X ℓ=1 ε(ω)−1 ε(ω) + ℓ+1 ℓ a2ℓ+1 R2ℓ+4 1 2ℓ(ℓ+ 1),(1.58b) with R the distance between the QE position rQE and the center of the spherical MNP. Moreover, Eind ( r = rQE, ω )can be obtained from Eq. (1.44) , and ˆαMNP ( ω ) ≡ α(ω)is given by Eq. (1.45). On the other hand, to obtain the optical response of a QE coupled to a spherical dimer, we numerically solve the quantities of Eq. (1.57) by following the methodology explained in Subsection 1.1.3. Thus, the optical response of the QE–MNPs system is obtained in this case in three steps: first, we calculate Eind ( r = rQE, ω )and ˆαMNP ( ω )by solving Eq. (1.35) for the external potential ϕext ( r, ω )given by Eq. (1.43) corresponding to plane-wave illumination. Second, ˆ G ( rQE,rQE, ω )and ˆαQE MNP ( ω )are obtained by solving the same equations but considering the external potential ϕext ( r, ω ) ≡ϕQE ( r )of a unitary point dipole given by Eq. (1.53) . Finally, we use Eq. (1.57) to calculate the total dipole moment induced at the coupled system, p(ω) = pMNP(ω) + pQE(ω). As an example, we show in Figure 1.7 the absorption cross-section spectra σabs ( ω )of an individual MNP (panel a) and a MNP dimer (panel b) interacting with a QE. The MNPs have a radius a = 5 nm and are characterized using the same Drude parameters as above (sodium). The QE is characterized [Eq. (1.45) ] 30 1.4. Plexcitonics: Quantum emitter excitons coupled to plasmons by an oscillator strength α0 QE = 2 au , an intrinsic loss rate γQE = 10 meV , and a transition frequency ωQE resonant with the main plasmon resonance of the nanostructures. In both situations, the presence of the QE strongly affects the optical absorption of the coupled system (blue lines) as compared to that of the isolated MNPs (dashed black lines). In particular, for the individual MNP Figure 1.7a shows a spectrally narrow Fano-like resonance at ω∼ωDP = 3 . 4 eV due to the destructive interference between the excitation of the plasmon and that of the exciton [ 173 , 174 ]. The effect of the QE on the optical response is stronger in the dimer configuration shown in Figure 1.7b due to the field enhancement in the gap and thus larger electromagnetic coupling between the QE exciton and the plasmonic resonances (see the following subsection for further details). In Chapters 3and 5, we study the impact of several quantum-mechanical phenomena in the optical response of QE–MNPs coupled systems. 1.4.3 Coupled harmonic-oscillator model The self-interaction Green’s function formalism explained in the previous subsection is generic and, in principle, can be applied to obtain the optical response of a QE interacting with any plasmonic nanostructure of arbitrary shape. However, using Eq. (1.57) requires the computation of the optical response of the plasmonic nanostructure (e.g. using numerical methods [ 112 – 115 ]), which can be difficult in many situations. Moreover, the exact numerical solution of the electromagnetic interaction between QEs and MNPs is usually difficult to interpret. In contrast, a simpler model based on coupled harmonic oscillators, explained in this subsection, can provide physical insight into such exciton–plasmon electromagnetic interaction [42,175–178]. In this model, the LSPP in the MNP and the exciton in the QE are each described as damped harmonic oscillators [ 42 ]. These oscillators represent, for example, the dipole moment induced at the MNP and at the QE, respectively, and are coupled through the electric near field induced by each structure. The equations of motion of the QE induced dipole moment, pQE ( t ), and that of the MNP, pMNP ( t ), are given under external plane-wave excitation Eext ( t ) = ReE0e−iωt in time domain by [42,178] ¨ pQE(t) + γQE ˙ pQE(t) + ω2 QEpQE(t) = α0 QEEext(t) + βMNPpMNP(t) ¨ pMNP(t) + γMNP ˙ pMNP(t) + ω2 MNPpMNP(t) = α0 MNPEext(t) + βQEpQE(t), (1.59) where ωMNP , γMNP and α0 MNP are the resonant frequency, intrinsic loss rate and oscillator strength associated with the single LSPP mode of the MNP, respectively. Only a single LSPP mode is considered in the MNP, so that a generalization of the model would be required to simultaneously account for the electromagnetic coupling between the exciton and multiple LSPPs, such as the BDP and BQP resonances supported by the metallic dimer shown in Figure 1.7. In Eq. (1.59) , the exciton–plasmon coupling is introduced via the electric near field induced by the 31 Chapter 1. Classical description of light–matter interaction MNP (QE), βMNPpMNP (βQEpQE), at the position of the QE (MNP). The direct excitation of the QE by the plane wave is usually much weaker than the excitation by the electric near field of the MNP, and thus one can consider α0 QEEext ( t ) = 0 in Eq. (1.59) [ 42 ]. Moreover, pMNP ( t )is expected to be much stronger than pQE ( t ). Under these assumptions, the expression of the approximated polarizability of the coupled system can be obtained in the frequency domain (using pMNP(t) = Re pMNP(ω)e−iωtand pQE(t) = Re pQE(ω)e−iωt) α(ω)≈α0 MNPω2 QE −ω2−iωγQE ω2 MNP −ω2−iωγMNPω2 QE −ω2−iωγQE−4g2ω2,(1.60) where we define the coupling strength g such that 4 g2ω2 = α0 QEα0 MNPβQEβMNP . The factor 4 is introduced so that g can be directly compared to the coupling strength used in cavity quantum-electrodynamics models [ 179 ]. This definition implies a dependence of g on ω , g∝ 1 /ω . However, in practice we consider g to be a constant (see below), consistent with the procedure adopted in the literature [ 42 ]. We expect that an alternative assumption for the definition of g would not substantially modify the results. The expression given by Eq. (1.60) can be useful e.g. to estimate the value of the coupling strength g between a QE exciton and a specific LSPP in the MNP. This can be done by fitting the exact polarizability α ( ω )of the QE–MNPs system rigorously obtained by using e.g. the methodology explained in Subsection 1.4.2 to the approximated result obtained from Eq. (1.60) . We follow this procedure in Chapter 5, which allows us to identify whether the studied QE–MNPs system is in the strong-coupling regime. In fact, the criteria adopted in the literature to identify the strong-coupling regime are intimately related to the value of g [ 180 ]. In brief, the less demanding criterion states that g > |γMNP −γQE|/ 4has to be fulfilled. Other criteria often used are g > ( γMNP + γQE ) / 4and the more restrictive g > (γMNP +γQE)/2. Figure 1.8 shows the absorption spectra σabs ( ω )[Eq. (1.49) ] of a QE–MNPs system obtained by solving Eq. (1.57) within the self-interaction Green’s function formalism (blue line) and the approximated value of σabs ( ω )obtained from Eq. (1.60) using the coupled harmonic-oscillator model (red dotted line). The parameters characterizing the system are the same as in Figure 1.7. In the case of the individual MNP (Figure 1.8a), the coupling strength g between the DP and the exciton is g = 18 meV < ( γMNP −γQE ) / 4 < ( γMNP + γQE ) / 4 < ( γMNP + γQE ) / 2, which indicates that the system is in the weak-coupling regime (here γMNP = γp ). For the metallic dimer (Figure 1.8b), the coupling strength g between the BDP and the exciton is g = 45 meV ≈ ( γMNP −γQE ) / 4 ≈ ( γMNP + γQE ) / 4, thus the system is nearly in the strong-coupling regime. For situations of strong coupling, the coupled harmonic-oscillator model depicted in this subsection can also be used to calculate the frequencies of the two polaritons resulting from the electromagnetic interaction between the exciton and the LSPP. Once the values of the coupling strength g is determined, as explained above, the 32 1.4. Plexcitonics: Quantum emitter excitons coupled to plasmons + Figure 1.8:(a) Absorption cross-section spectra σabs ( ω )of a QE–MNP coupled system obtained by using the self-interaction Green’s function formalism explained in Subsection 1.4.2 (blue line) and a coupled harmonic-oscillator model (Eq. (1.60) , red dotted line). The MNP has a radius a = 5 nm , and it is characterized by a Drude dielectric function with ωp = 5 . 89 eV and γp = 0 . 21 eV . Within the coupled-oscillator model ωMNP = ωDP = 3 . 4 eV is used. The QE, located at 1 . 5 nm from the surface of the MNP, is characterized by an oscillator strength α0 QE = 2 au , an intrinsic loss rate γQE = 10 meV , and a transition frequency ωQE = ωDP in resonance with the DP of the MNP. The coupling strength is g = 18 meV , and it is obtained by fitting the exact result obtained from Eq. (1.57) to the expression given by Eq. (1.60) .(b) Same as in panel (a) but for a MNP dimer characterized by a gap separation D = 3 nm . Within the coupled-oscillator model ωMNP = ωBDP = 3 . 15 eV is used. The QE emitter is at resonance with the BDP mode of the dimer, ωQE = ωBDP , and the coupling strength is g = 45 meV . In both panels, γMNP =γpis used within the coupled-oscillator model. frequencies of the upper ( ω+ ) and the lower ( ω− ) polaritons can be obtained by solving the eigenmodes of Eq. (1.59) [179] ω±=1 2(ωMNP +ωQE)±1 2Re  s4g2+ωMNP −ωQE +iγQE −γMNP 22   . (1.61) According to Eqs. (1.61) , in a resonant exciton–plasmon system ( ωMNP ≈ωQE ), ω+ and ω− are separated by a factor 2 qg2−(γQE −γMNP)2/16 as long as g > |γQE −γMNP|/ 4, which justifies the less demanding criterion of strong coupling mentioned above. Moreover, for situations where g >> |γQE −γMNP|/ 4, the separation is 2 g , leading to the so-called Rabi splitting between the upper and lower polaritons in strongly coupled systems [ 181 – 183 ]. In contrast, in the weakcoupling regime, where g < |γQE −γMNP|/ 4, there is no splitting between ω+ and ω− , indicating that no polaritonic (hybrid) modes are created as a consequence of the electromagnetic interaction between the exciton and the LSPP. 33 Chapter 1. Classical description of light–matter interaction 1.5 Summary In summary, we have presented in this chapter the classical viewpoint of plasmonic excitations supported by MNPs and their electromagnetic coupling to QE excitons. First, we have reviewed Maxwell’s equations in the linear-response regime, and focused on the local-response approximation (LRA) that neglects the spatial nonlocality of the metal response. Next, we have introduced the nonretarded approximation adopted in this thesis and presented a method to obtain the optical response of spherical individual nanoparticles and dimers. Then, after briefly describing bulk and surface plasmons, we have paid special attention to the excitation of LSPPs in plasmonic nanostructures of finite geometry such as spherical MNPs and dimers. We have also introduced a simple model to understand the nonlinear optical response from MNPs. Finally, we have presented the theoretical approaches used in this thesis to describe the electromagnetic interaction between QEs and MNPs within a classical linear-response framework. The concepts explained in this chapter thus provide the grounds of classical light– matter interaction and serve as a reference to study quantum phenomena that are out of the reach of classical descriptions. In this thesis, quantum effects arising in the light–matter interaction are studied within TDDFT as explained in the following Chapter 2. 34 Chapter 2 QUANTUM MANY-BODY DESCRIPTION OF LIGHT–MATTER INTERACTION The classical theoretical framework introduced in Chapter 1can be used to accurately describe light–matter interaction when the characteristic dimensions of the system, such as the size of the MNP or the gap separation in MNP ensembles, are relatively large. However, when small MNPs ( ≲ 10 nm ) or ultra-narrow gaps ( ≲ 1 nm ) are considered, the quantum nature of the electrons dynamics becomes important and classical descriptions are no longer valid. In this context, a methodology capable of describing the electronic structure from a quantum manybody perspective is required, which is a considerable challenge [ 184 ]. For example, to describe the ground state of a small sodium MNP containing 1000 atoms (with 11 electrons per atom), one would have to solve the Schrödinger equation for a wave function depending on 33000 spatial variables (three spatial variables per electron without considering spin degeneracy and neglecting the degrees of freedom of the nuclei). Solving this tremendously complex problem is out of current computational capabilities [185]. To reduce the computational complexity of this many-body problem, one can use density-functional theory (DFT), a rigorous formalism that deals with the ground-state electron density rather than with the many-electron wave function. DFT can be applied only to study ground-state electronic properties, so that to determine electronic excitations in metals (and other materials) resulting from light– matter interaction it is necessary to adopt the time-dependent extension of DFT, the so-called time-dependent density functional theory (TDDFT). This approach addresses the time evolution of the electron density when the system is subjected to a time-dependent external potential. On the other hand, a semiclassical approach 35 Chapter 2. Quantum many-body description of light–matter interaction referred to as the surface-response formalism (SRF) first introduced by Peter Feibelman in the 1980s [ 98 ] which incorporates the parameters d⊥ ( ω )and d∥ ( ω ), has prompted great interest and practical use in the Nanophotonics community over the last few years to account for certain quantum many-body effects in the optical response of metals. The advantage of using the SRF is that it is much simpler computationally as compared to TDDFT, and thus allows for studying quantum effects in larger plasmonic nanostructures. However, the SRF is a cruder approximation than TDDFT, and cannot account for all quantum many-body phenomena. In Sections 2.1 and 2.2, we briefly recall the fundamentals of DFT and TDDFT, and describe the corresponding quantum many-body algorithms based on the wave-packet propagation (WPP) method employed in this thesis to address the dynamics of the electron density in time domain. We apply this TDDFT framework in Sections 2.3 and 2.4 to introduce some of the main quantum many-body phenomena manifested in the linear and nonlinear optical response of canonical plasmonic systems such as individual spherical and cylindrical MNPs and their dimers. Finally, we explain in Section 2.5 the fundamentals of the SRF used in this thesis (complementary to TDDFT) to provide additional insights on the study of particular quantum surface effects on the optical response. 2.1 Fundamentals of density functional theory (DFT) Density-functional theory (DFT) allows us to determine the ground-state electronic properties of an interacting many-electron system by only calculating the groundstate (or equilibrium) electron density n0 ( r ), i.e., without the need of the exact wave function satisfying the many-electron Schrödinger equation. According to the work by Hohenberg and Kohn [ 186 ], n0 ( r )completely determines the ground-state energy, E0 , and all other electronic properties of the many-electron system subjected to an external time-independent potential. This potential can be, for example, the attractive Coulomb potential Vion ( r )created by the positively charged ions in a metal, as we consider here. In this case, the ground-state energy E0 of a metal can be expressed as a unique functional of n0(r)(denoted by the square brackets), E0[n0(r)] = Zdrn0(r)Vion(r) + 1 2Z Z drdr′n0(r)n0(r′) |r−r′|+G[n0(r)],(2.1) where the first and second terms on the right-hand side (RHS) are the energy due to electron–ion and electron–electron Coulomb interaction in a metal, and G [ n0 ( r )] is a universal functional of the density n0 ( r )valid for any number of electrons and any potential Vion ( r ). G [ n0 ( r )] accounts for the kinetic energy and exchange–correlation energy (associated e.g. with the Pauli exclusion principle) of the interacting many-electron system. If G [ n0 ( r )] were known, determining the ground-state energy E0 and electron density n0 ( r )of a many-electron system could 36 2.1. Fundamentals of density functional theory (DFT) be directly obtained by minimizing E0 [ n0 ( r )] relative to n0 ( r )according to the variational principle [ 186 ]. Thus, the work by Hohenberg and Kohn completely changes the paradigm of the electronic many-body problem, since dealing with the electron density n0 ( r )as the fundamental quantity instead of the the many-electron wave function allows the study of electronic properties of many-electron systems intractable in the past. Notice that Eq. (2.1) is formally exact. Unfortunately, the universal functional G [ n0 ( r )] is generally unknown and, therefore, DFT becomes in practice an approximation. Considerable effort has been devoted over the years to find suitable approximations for G[n0(r)] [187–189]. Kohn and Sham [ 190 ] proposed to express G [ n0 ( r )] as a sum of two functionals, G[n0(r)] = Ts[n0(r)] + Exc[n0(r)],(2.2) where Ts[n0(r)] = X j∈occ ZdrΨ0 j(r)∗ˆ TΨ0 j(r)(2.3) is the kinetic energy of an auxiliary system of non-interacting electrons 6 and Exc [ n0 ( r )] is the exchange–correlation energy that contains all remaining manybody interactions. In Eq. (2.3) , ˆ T = −1 2∇2 is the kinetic-energy operator, and the summation runs over the occupied ( j∈occ ) time-independent Kohn-Sham (KS) orbitals Ψ 0 j ( r )that determine the equilibrium electron density n0 ( r )of the many-body system, n0(r) = X j∈occ χj|Ψ0 j(r)|2,(2.4) with the statistical factors χj accounting for both spin and symmetry degeneracy. Note that, according to the definition of Eq. (2.4) , n0 ( r )is considered to be positive. The application of the variational principle to Eq. (2.1) using the functionals given by Eqs. (2.2) and (2.3) results in the following equation [191]: ˆ T+Vion(r) + Zdr′n0(r′) |r−r′|+δExc[n0(r)] δn0(r) | {z } ˆ H0[n0(r)] Ψ0 j(r) = ϵjΨ0 j(r),(2.5) where Ψ0 j(r)and ϵjare the time-independent KS orbitals and energies7. Thus, the equilibrium electronic density n0 ( r )[Eq. (2.4) ] of the true many-body system can be obtained within the KS scheme from the solutions of the Schrödinger equation [Eq. (2.5) ] of an auxiliary system of non-interacting electrons using an 6Ts is thus an explicit functional of the Kohn–Sham orbitals Ψ 0 j ( r ), but an implicit functional of n0(r)according to Eq. (2.4) [72]. 7 Note that Ψ 0 j ( r )and ϵj in Eq. (2.5) are the eigenfunctions and eigenvalues of the auxiliary non-interacting electron system used to construct the exact electron density n0 ( r ), and therefore, ϵj and Ψ 0 j ( r )have no direct physical interpretation [ 192 , 193 ]. An exception is made for the energy of the highest occupied KS orbital, which can be used e.g. to estimate the work function of a metal [194]. 37 Chapter 2. Quantum many-body description of light–matter interaction effective time-independent Hamiltonian ˆ H0[n0(r)] = ˆ T+Veff[n0(r)],(2.6) where Veff[n0]is the effective one-electron potential given by Veff[n0(r)] = Vion(r) + Zdr′n0(r′) |r−r′|+Vxc[n0(r)].(2.7) In Eq. (2.7) , the second term on the RHS is the electron–electron potential due to Coulomb interaction in the single-electron picture, and Vxc [ n0 ( r )] is the exchange– correlation potential obtained from Exc[n0(r)] through the relationship Vxc[n0(r)] = δExc[n0(r)] δn0(r),(2.8) which accounts for all many-body interactions that are not present in ˆ T and in the second term of Eq. (2.7) , but exist in the true interacting many-electron system. According to the previous discussion about G [ n0 ( r )], the exchange–correlation energy functional Exc [ n0 ]is not known exactly. We explain in Subsection 2.1.1 the approximation adopted in this thesis to calculate Exc[n0(r)]. Equations (2.4) and (2.5) are referred to as the time-independent KS equations, and have to be solved self-consistently: one can start with an assumed electron density n0 ( r ), then construct Veff[n0(r)] from Eq. (2.7) , and finally obtain a new value of n0 ( r )using Eqs. (2.5) and (2.4) . This procedure is repeated until convergence of the results is achieved. In the following Subsections 2.1.1 and 2.1.2, we describe the approximations used in this thesis to compute the potentials Vxc[n0(r)] and Vion(r). 2.1.1 The local-density approximation (LDA) The simplest way to determine the exchange–correlation energy functional Exc [ n0 ( r )] in Eq. (2.8) is to adopt the local-density approximation (LDA). If we assume that n0 ( r )does not change rapidly, the variation of Exc [ n0 ( r )] with respect to the gradient of n0(r)can be neglected, and Exc[n0(r)] can be expressed as [62] Exc[n0(r)] = Zdrn0(r)ϵxc(n=n0(r)),(2.9) where ϵxc ( n )is the exchange and correlation energy per electron of a homogeneous electron gas with average electron density n [ 190 ]. Thus, within the LDA, the exchange–correlation energy per particle located at position r in an inhomogeneous system with density n0 ( r )is approximated by the exchange–correlation energy per particle of a uniform electron gas with the same density, n = n0 ( r ). The LDA has been widely used to determine the electronic properties of many systems, including atoms and solids where the density does not vary slowly. For metallic surfaces, 38 2.1. Fundamentals of density functional theory (DFT) the LDA is found to provide accurate results, which might be surprising since the equilibrium density n0(r)varies rapidly near the metal surface [62]. Using Eq. (2.8) , we obtain the exchange–correlation potential within the LDA 8 , Vxc[n0(r)] = ∂n0(r)ϵxc(n=n0(r)) ∂n0(r)=ϵxc(n=n0(r)) + n0(r)∂ϵxc(¯n) ∂n n=n0(r) , (2.10) which requires an analytical expression of ϵxc ( n ). Several approximations have been proposed to that end [ 191 , 195 – 197 ], and throughout this thesis, we use the exchange–correlation energy-density functional ϵxc ( n )given by Gunnarsson and Lundquist [198]9. ϵxc(n) = −1 20.916/rs+ 0.0666(1 + x3)ln(1 + 1/x)−x2+x 2−1/3),(2.11) with rs = 3 4π¯n1/3 the Wigner–Seitz radius [Eq. (1.11) ], and x = rs/ 11 . 9. Finally, according to Eq. (2.10) , the exchange–correlation potential that we use in this thesis is given by [199] Vxc[n0(r)] = −1 21.222/rs+ 0.0666 ln(1 + 11.4/rs)n=n0(r) .(2.12) 2.1.2 The jellium model of free-electron metals To obtain the potential Vion ( r )in Eq. (2.7) , we adopt in this thesis the jellium model of free-electron metals [ 194 , 200 , 201 ], where the ions at the lattice sites are modeled as a uniform positive background charge with density n+(r) =    ninside the metal 0outside the metal .(2.13) Thus, within the jellium model, n (or, equivalently, rs [Eq. (1.11) ]) is the only parameter needed to characterize the metal. For example, Al is modeled with rs= 2.07 a0, Na with rs= 4 a0, and K with rs= 4.96 a0[62]. This positive charge density n+(r)creates an attractive Coulomb potential Vion(r) = −Zdr′n+(r′) |r−r′|.(2.14) 8 It can be shown that, for a functional F [ ρ ( r )] of the form F [ ρ ( r )] = Rdrfr, ρ ( r )  , the functional derivative is given by δF δρ(r) = ∂f ∂ρ . A demonstration can be found in https://www.youtube.com/watch?v=_ntUQ_WBp0U. 9 We consider the spin-unpolarized case with ξ = 0. Moreover, note that in ref. 198 the expression for ϵxc ( n )is given in Rydberg atomic units (while we use Hartree atomic units) and thus it differs from our Eq. (2.11) by a factor 1 2. 39 Chapter 2. Quantum many-body description of light–matter interaction from the time-to-frequency Fourier transform F[Eq. (1.3)], p(ω) = ZTf 0 dt p(t)eiωt F(t),(2.36a) Eind(r, ω) = ZTf 0 dt Eind(r, t)eiωt F(t),(2.36b) j(r, ω) = ZTf 0 dt j(r, t)eiωt F(t),(2.36c) where Tf is the total propagation time used in our simulations, which must be long enough to achieve convergence. Importantly, the time-dependent function F ( t )(or filter)inEq. (2.36) is introduced to attenuate the collective time-dependent charge-density oscillations (losses), since the ALDA-TDDFT scheme adopted in this thesis does not account for decay and dephasing processes related to the scattering of electrons with phonons, nor for the intrinsic losses due to inelastic electron–electron interactions. This is a well-known failure of the ALDA used for the exchange–correlation potentials Vxc [ n ( r, t )] [ 72 , 220 , 241 , 242 ], and thus a filter F ( t )has to be applied to mimic such intrinsic losses in the system. In practice, in this thesis we employ two different filters F ( t )depending on the external excitation, as described in Sections 2.3 and 2.4. Calculation of the projected density of electronic states (PDOS) within the WPP method Given an equilibrium electron density n0 ( r )and a ground-state effective potential Veff ( r )of a system [Eq. (2.17) ], it is possible to access the energies of both the occupied and unoccupied KS one-electron states of the system by analyzing the projected density of electronic states (PDOS), Σ( ω ). Importantly, the PDOS can be obtained within the WPP method described here by propagating an initial wave packet Φ( r, t = 0) = Φ 0 ( r )of a particular symmetry according to the time-dependent Schrödinger equation under the (time-independent) one-electron Hamiltonian ˆ H = ˆ T + Veff ( r ). This Hamiltonian corresponds to the Hamiltonian of a single electron subjected to the effective ground-state potential Veff ( r )of the nanostructure under study. In this case, using the WPP method to obtain the PDOS we do not apply any external potential Vext ( r, t ), and Veff ( r )does not vary in time. However, the PDOS provides information about the occupied and unoccupied electronic states that will be involved in electronic transitions excited optically when an external time-dependent potential is applied. The PDOS Σ(ω)projected onto the initial wave packet Φ0(r)is given by Σ(ω) = ∞ X j=1 |cj|2δ(ω−ϵj),(2.37) with δ ( ω−ϵj )the Dirac delta, ϵj the eigenenergies of the one-electron Hamiltonian 46 2.2. Fundamentals of time-dependent density functional theory (TDDFT) ˆ H=ˆ T+Veff(r), and cj=⟨ϕj(r)|Φ0(r)⟩(2.38) the complex coefficients corresponding to the expansion of the initial wave packet Φ0(r)into the eigenfunctions ϕj(r)of the Hamiltonian ˆ H, Φ0(r) = ∞ X j=1 cjϕj(r).(2.39) The eigenfunctions ϕj ( r )are not known a priori, and we apply the WPP method to obtain Σ(ω)[Eq. (2.37)] without the need to calculate ϕj(r). The PDOS Σ( ω )given by Eq. (2.37) represents the number of one-electron states of a particular spatial symmetry that the nanostructure sustains at a given energy level ω and spatial region. The PDOS includes the contribution of both the occupied and the unoccupied states, and it gives information about the degree of localization of a particular electronic state at a certain spatial region (determined by the initial wave packet Φ 0 ( r )used in the propagation). Notice that, since the PDOS shows the one-electron energy states available in a (fictitious) non-interacting KS electron system (see Section 2.1), the resonant energies optically excited at the interacting many-electron system will be renormalized with respect to the non-interacting one via Hartree and exchange–correlation potentials [ 192 , 193 ]. However, the PDOS still provides useful insights into the properties of the electronic structure of the system. For example, we study in Chapter 5the PDOS in the electronically coupled QE–MNPs system to quantify the degree of electronic hybridization of the occupied and unoccupied electronic states. In order to calculate Σ( ω )[ 229 , 243 ], we choose an initial wave packet Φ( r, t =0)=Φ 0 ( r )of a particular symmetry and propagate it according to the time-dependent Schrödinger equation corresponding to ˆ H=ˆ T+Veff(r), Φ(r, t) = e−i(ˆ T+Veff(r))Φ0(r),(2.40) by using the WPP algorithm described above [Eqs. (2.31) and (2.32) ]. Note that no external perturbation is used in the WPP calculations of the PDOS, Vext ( r, t ) = 0, and that the self-consistent procedure is not applied so that Veff ( r )remains constant in time. Once the time evolution of the initial wave packet Φ( r, t )is obtained by solving Eq. (2.40) , to calculate Σ( ω )we first apply the time-to-frequency Laplace transform ˆ Lωto Φ(r, t), ˆ LωΦ(r, t) = Z∞ 0 dt ei(ω+iζ)tΦ(r, t) = Z∞ 0 dt ei(ω−ˆ H+iζ)tΦ0(r) = i (ω−ˆ H+iζ)Φ0(r), (2.41) where ζ→ 0 + is a small positive number. Using Eq. (2.39) ,Eq. (2.41) can be 47 Chapter 2. Quantum many-body description of light–matter interaction written as ˆ LωΦ(r, t) = i∞ X j=1 cjϕj(r) (ω−ϵj+iζ).(2.42) Then, projecting Eq. (2.42) onto the initial wave packet Φ0(r)one obtains ⟨Φ0(r)|ˆ LωΦ(r, t)⟩=i∞ X j=1 |cj|2 (ω−ϵj+iζ).(2.43) Finally, applying the Sokhotski–Plemelj theorem to Eq. (2.43) , the PDOS Σ( ω ) given by Eq. (2.37) can be expressed as Σ(ω) = 1 πlim ζ→0+Re{⟨Φ0(r)|ˆ LωΦ(r, t)⟩}, =1 πlim ζ→0+Re{ˆ Lω⟨Φ0(r)|Φ(r, t)⟩} | {z } A(t) ,(2.44) where A ( t ) = ⟨ Φ 0 ( r ) | Φ( r, t ) ⟩ is the autocorrelation function that can be calculated in time domain using the WPP method. Notice that cj are the coefficients of Φ 0 ( r )decomposed into the eigenfunctions ϕj ( r )of the Hamiltonian ˆ H [Eq. (2.38) ], so that choosing different initial wave packets Φ 0 ( r )leads to a different PDOS Σ( ω )[ 244 ]. This can be useful to focus on specific electronic states with a special symmetry or spatial distribution (as we do in Chapter 5), since only these states of the system that are not orthogonal to the initial wave packet Φ 0 ( r )can be accessed by the WPP method used to calculate the PDOS. 2.3 Linear optical response of canonical plasmonic nanostructures addressed within TDDFT In this section, we analyze TDDFT results of canonical systems to illustrate the general properties of the linear optical response of plasmonic nanostructures. In particular, we discuss in Subsection 2.3.1 the influence of the nanoparticle size in the optical response of spherical MNPs described within the jellium model, and introduce the quantum-mechanical concepts of Friedel oscillations, electron spill out, and surface-enabled Landau damping. In Subsection 2.3.2, we analyze the effect of the size and the gap separation in the optoelectronic response of spherical MNP dimers, and introduce the concept of electron tunneling. Finally, we study in Subsection 2.3.3 the linear optical response of individual cylindrical metallic nanowires and their dimers, which allows us to study larger nanostructures than the spherical MNPs because the symmetry of the system reduces the computational demands. 48 2.3. Linear optical response of canonical plasmonic nanostructures addressed within TDDFT 2.3.1 Individual spherical metallic nanoparticles In this subsection, we consider individual spherical MNPs surrounded by vacuum. The electronic structure of the system is described within the jellium model of free-electron metals (see Subsection 2.1.2) using a Wigner–Seitz radius rs = 4 a0 that corresponds to sodium. We consider closed-shell MNPs, and the number of conduction electrons Ne is varied to study the influence of the MNP size in the optical response. The radius of the MNPs, a, is determined from a=N1/3 ers.(2.45) Before showing the results, we provide the main numerical details to carry out the TDDFT simulations. Numerical implementation To obtain the ground-state KS orbitals of individual spherical MNPs containing Ne conduction electrons using DFT, we take advantage of the spherical symmetry of the problem and solve the time-independent KS equations given by Eqs. (2.4) and (2.5) in a spherical coordinate system r = {r, θ, φ} (Figure 1.2a). We write the time-independent KS orbitals as Ψ0 jr,ℓ,m(r)≡Ψ0 jr,ℓ,m(r, θ, φ) = 1 rψ0 jr,ℓ(r)Ym ℓ(θ, φ),(2.46) where the sequence of the radial ( jr = 1 , 2 , . . . )and angular ( ℓ = 0 , 1 , . . . )quantum numbers is limited by the condition ϵjr,ℓ ≤EF , where ϵjr,ℓ is the energy of the ground-state KS orbital (see below) and EF is the Fermi energy. The magnetic quantum number can take the values m = ( −ℓ, . . . , 0 , . . . , ℓ ). In Eq. (2.46) , the radial part ψ0 jr,ℓ ( r )satisfies the following one-dimensional Schrödinger-like equation [see Eqs. (2.5) and (1.24)]: −1 2 d2 dr2+ℓ(ℓ+ 1) 2r2+Veff[n0(r)] | {z } ˆ H0[n0(r)] ψ0 jr,ℓ(r) = ϵjr,ℓψ0 jr,ℓ(r),(2.47) with Veff [ n0 ( r )] given by Eq. (2.17) , and the equilibrium electron density n0 ( r )is given by n0(r)=2 X jr,ℓ∈occ 1 4π 2ℓ+ 1 r2|ψ0 jr,ℓ(r)|2.(2.48) In Eq. (2.48) , the factor 2is due to the spin, and (2ℓ+ 1) due to the degeneracy of the KS orbital with the same angular quantum number, ℓ , and different magnetic quantum number, m . Because of the spherical symmetry of the problem, n0 ( r ) and Veff[n0(r)] depend only on the radial coordinate r. The Hartree potential VH [ n0 ( r )] contained in Veff [ n0 ( r )] [Eq. (2.17) ] is calculated from Poisson’s equation [Eq. (2.16) ] by defining ˜ VH [ n0 ( r )] = rVH [ n0 ( r )]. This 49 Chapter 2. Quantum many-body description of light–matter interaction definition allows us to obtain the Hartree potential from d2 dr2˜ VH[n0(r)] = −r4π(n0(r)−n+(r)),(2.49) which can be solved by directly applying the inverse matrix of operator d2 dr2 to the RHS of Eq. (2.49) using space-to-momentum Fourier sine transform. As described in Section 2.1, the time-independent KS equations given by Eqs. (2.47) and (2.48) are solved self-consistently using an iterative procedure by diagonalization of the Fourier grid Hamiltonian (FGH) [ 237 ] obtained from H0 [ n0 ( r )] [Eq. (2.47) ], where the KS orbitals ψ0 jr,ℓ ( r )are represented in a real-space mesh of equidistant points in the coordinate r , and the space-to-momentum Fourier sine transform is used to compute the operator d2 dr2 [ 230 , 232 ]. At the first iteration, the effective potential is initialized to Veff(r) = V0 eff e−(r−a−1) 1 + e−(r−a−1) ,(2.50) where V0 eff is a (negative-valued) parameter that determines the depth of the potential. More details on the self-consistent procedure can be found in ref. 245. Once the ground-state KS orbitals are calculated using DFT, we obtain the linear optical response of the spherical MNP within TDDFT by applying the WPP method described in Subsection 2.2.1. We consider here the following timedependent external potential, Vext(r, t) = δ(t)E0∆t r cos θ, (2.51) which represents a perturbation at the initial time t = 0 within the dipole approximation corresponding to a plane-wave electric field polarized along the z -axis. In Eq. (2.51) , δ ( t )is the Dirac delta function, E0 is the amplitude of the external perturbation (we typically use E0∼ 10 −5au , weak enough so that the linear-response approximation holds), and ∆ t is the propagation time step used in our simulations, typically ∆ t∼ 0 . 25 − 0 . 1 au (1 au ≈ 2 . 419 × 10 −2fs , see Appendix A). The external potential given by Eq. (2.51) breaks the spherical symmetry of the system but preserves the rotational symmetry with respect to the z -axis (azimuthal symmetry), and therefore the magnetic quantum number m is still a good quantum number. It is thus convenient to express the time-dependent KS orbitals as15 Ψj,m(r, θ, φ, t) = 1 rψj,m(r, θ, t)1 √2πeimφ,(2.52) where, here, the quantum number j replaces the pair {jr, ℓ} defining the ground15 In this thesis, all the considered external potentials Vext ( r, t )acting on spherical MNPs preserve the rotational symmetry of the system and thus we can always write the KS orbitals of spherical MNPs using Eq. (2.52) and follow the procedure described in this subsection. 50 2.3. Linear optical response of canonical plasmonic nanostructures addressed within TDDFT state KS orbital. The orbitals ψj,m ( r, θ, t )in Eq. (2.52) are obtained from the time-dependent KS equations [Eq. (2.19) ], with the kinetic-energy operator ˆ T taking the form [Eq. (1.24)] ˆ T=−1 2∂2 ∂r2+1 r21 sinθ ∂ ∂θ sin θ∂ ∂θ−m2 sin2θ.(2.53) The time evolution of ψj,m ( r, θ, t )and ψj,−m ( r, θ, t )is identical, and therefore we only propagate ψj,m ( r, θ, t )for m≥ 0. We can thus express the time-dependent electron density as n(r, θ, t)X {j,m≥0}∈occ 1 2πχm 1 r2|ψj,m(r, θ, t)|2,(2.54) with χm=   2for m= 0 4for m > 0 (2.55) accounting for both the spin and ±m degeneracy. The orbitals ψj,m ( r, θ, t )are propagated in time by adapting the WPP scheme given by Eq. (2.32) to the use of spherical coordinates. We represent the KS orbitals ψj,m ( r, θ, t )in a meshgrid extending in radial direction up to 25 − 35 a0 from the jellium edge, using a constant radial spacing of ∆ r∼ 0 . 35 − 0 . 5 a0 . The angular variable θ is discretized from 0 to πusing 60 −120 points. In particular, the operator e−i∆t−1 2 ∂2 ∂r2 linked to the first term on the RHS of Eq. (2.53) is applied to ψj,m ( r, θ, t )in the reciprocal space using space-to-momentum Fourier sine transforms. On the other hand, the operator e−i∆t−1 2r21 sin θ ∂ ∂θ (sin θ∂ ∂θ )−m2 sin2θ linked to the second and third terms on the RHS of Eq. (2.53) is applied in real space by expanding ψj,m ( r, θ, t )in a basis of associated Legendre polynomials Pm ℓ ( cosθ ), taking advantage of the fact that Pm ℓ ( cosθ )are the eigenfunctions of the operator n1 sin θ ∂ ∂θ sin θ∂ ∂θ −m2 sin2θo . This procedure is described in detail in ref. 75. The linear optical response of the MNP is analyzed by calculating its absorption cross-section spectrum, σabs ( ω ) = 4π cIm{α ( ω ) } [Eq. (1.49) ], with α ( ω )the polarizability of the system calculated from α(ω) = 1 E0∆tp(ω)16. p(ω)is obtained using the time-to-frequency Fourier transform given by Eq. (2.36a) by considering 16 p ( ω )is the frequency-dependent dipole moment induced along the z -axis, obtained according to Eqs. (2.33) and (2.36a). 51 Chapter 2. Quantum many-body description of light–matter interaction 2.5 3 3.5 4 Figure 2.1:(a) Time evolution of the dipole moment p ( t )induced at a spherical sodium MNP ( rs = 4 a0 ) containing Ne = 1074 electrons in response to the external excitation given by Eq. (2.51) . The results directly obtained from TDDFT within the WPP method (gray line) are damped by using different values of the attenuation factor η [Eq. (2.56) ], as indicated in the legend. (b) Absorption cross-section spectra σabs ( ω )of the same MNP as in panel (a) obtained for different values of the attenuation factor η . Using different values of the factor η (with η < κ ) leads to very similar widths of the DP resonance, κ . The results corresponding to η = 0 are divided by 2 for clarity. in this case the filter F(t) = e−η/2t,(2.56) with an attenuation factor η∼ 0 . 05 − 0 . 2 eV smaller than the typical plasmon resonance width (as determined e.g. through the classical calculations in Chapter 1,γpin Eq. (1.10)). We show in Figure 2.1a the effect of the filter F ( t )given by Eq. (2.56) in the result of the time-dependent dipole moment p ( t )induced at a spherical MNP containing Ne = 1074 conduction electrons. Without applying the filter ( η = 0, gray line), p ( t )oscillates in time showing revivals because, as discussed in Subsection 2.2.1, the ALDA-TDDFT model adopted in this thesis does not include intrinsic dissipation processes in the free-electron gas. When an attenuation factor η > 0 is applied, p ( t )is exponentially damped, thus accounting for the losses of the metal phenomenologically. This attenuation factor η also affects dramatically the absorption cross-section spectrum σabs ( ω )of the system as illustrated in Figure 2.1b. Indeed, for η = 0 the plasmon peak in σabs ( ω )at around ω = 3 . 2 eV is fragmented into a set of discrete lines associated with the single electron–hole excitations that build up the plasmon 17 (see below) [ 206 , 246 ]. The inclusion of η broadens these single electron–hole excitation peaks and allows us to retrieve the broad Lorentzian-like profile of the plasmon resonance similarly to the classical 17 The width of the peaks associated to single electron–hole excitations in Figure 2.1b for η = 0 is due to the finite calculation time Tf , where in this case we use Tf∼ 13000 au . However, we typically use Tf∼ 3000 − 4000 au in this thesis for linear-response calculations, enough to achieve convergence when using η= 0.05 −0.2eV. 52 2.3. Linear optical response of canonical plasmonic nanostructures addressed within TDDFT prediction shown in Figure 1.4d. We employ the filter given by Eq. (2.56) in Chapter 3, Chapter 4, and Chapter 5, where we study the linear optical response of different plasmonic systems. Figure 2.2:(a) Equilibrium electron density n0 ( r )(solid lines) and background jellium denisty n+ ( r )(dashed filled lines) as a function of the radial coordinate, r , for spherical MNPs characterized by a Wigner–Seitz radius rs = 4 a0 (sodium) and containing different numbers of conduction electrons: Ne = 338 electrons (purple), 638 electrons (brown), 1074 electrons (blue), 2260 electrons (green), and 4458 electrons (red). Results are normalized to the average electron density ¯n = 4 3πr3 s−1 [Eq. (1.11) ]. (b) Induced electron density δn ( r, ω )(multiplied by r2 ) along the radial axis r ( θ = 0) at the DP frequency ω = ωDP = 3 . 25 eV for a sodium MNP containing 4458 electrons in response to the external plane-wave excitation given by Eq. (2.51) . The dashed red line represents the position of the jellium edge of the MNP. An attenuation factor η= 0.2eV [Eq. (2.56)] is used to perform the Fourier transform. Influence of the size of small spherical nanoparticles on their optical response We show in Figure 2.2a the equilibrium electron density profile n0 ( r )obtained for spherical MNPs containing different number of conduction electrons within the range Ne = 338 − 4458, resulting in a radius a∼ 1 . 5 − 3 . 5 nm [Eq. (2.45) ]. Quantummechanical phenomena such as electron spill-out [ 247 ] and Friedel oscillations [ 248 , 249 ] of the equilibrium electron density n0 ( r )are observed. The electron spill out of n0 ( r )is a consequence of the finite potential barrier at the MNP surface that allows electrons to spread outside the background jellium edge located at r = a . Moreover, Friedel oscillations are due to electron reflection at r = 0 and at the MNP boundary, and can be understood as a manifestation of the Gibbs phenomenon occurring in the Fourier series of step-like functions: the electron density n0 ( r )is given by a summation of a finite number of ψ0 jr,ℓ ( r )orbitals [Eq. (2.48) ] and thus exhibits oscillations along the radial coordinate r . In addition, n0 ( r )in Figure 2.2a features either a peak or a dip at r = 0 depending on the number of conduction electrons. While for Ne = 338 electrons (purple), 638 electrons (brown) and 1074 53 Chapter 2. Quantum many-body description of light–matter interaction electrons (blue) there is a dip at r = 0, for Ne = 2260 electrons (green) and 4458 electrons (red) there is a peak. This behavior is related to the relative contribution of the ψ0 jr,ℓ ( r )orbitals with ℓ = 0 (which depends upon the number of occupied closed shells), since only those with ℓ = 0 can contribute to the electron density n0 ( r )precisely at r = 0 and only lowℓ orbitals ψ0 jr,ℓ ( r )contribute close to the center of the MNP because of the centrifugal potential ℓ ( ℓ + 1) /r2 [Eq. (2.47) ] that forbids electrons with high angular quantum number ℓ to approach the center of the MNP. Figure 2.2b shows the electron density δn ( r, ω )induced at the dipolar plasmon (DP) frequency ωDP = 3 . 25 eV in a spherical MNP ( Ne = 4458 is considered as an example) in response to an external plane-wave excitation. The induced density δn ( r, ω )is predominantly located near the surface of the MNP as expected from classical descriptions. However, due to the nonlocal dynamic screening of conduction electrons described within the jellium model [ 58 , 126 ], δn ( r, ω )spreads beyond the limits of the classical sharp edge at r = a (not to be confused with the spill out of the equilibrium electron density n0(r)shown in Figure 2.2a). 0 100 200 300 400 500 2.5 3 3.5 4 2.5 3 3.5 4 Figure 2.3:(a) TDDFT results of the absorption cross-section spectra σabs ( ω )of spherical MNPs with different number of conduction electrons, Ne =338 (purple), 638 (brown), 1074 (blue), 2260 (green), and 4458 (red) . The electronic structure of the MNPs is described within the jellium model using a Wigner–Seitz radius rs = 4 a0 that corresponds to sodium. An attenuation parameter η = 0 . 07 eV is used [Eq. (2.56) ]. (b) Classical LRA results of σabs ( ω )for spherical MNPs with different radius a , as obtained from Eq. (1.45) . a = 27 . 9 a0(purple), 34 . 4 a0(brown), 40 . 96 a0(blue), 52 . 49 a0(green), and 65 . 83 a0(red) . The value of the radius a is determined according to Eq. (2.45) . A Drude dielectric function [Eq. (1.10) ] with ωp = 5 . 63 eV and γp = 0 . 175 eV is used for all MNP sizes, which reproduces the TDDFT data for the largest MNP (red curve). We next show in Figure 2.3 the absorption cross-section spectra σabs ( ω )of the spherical sodium MNPs. The TDDFT results (panel a) illustrate that quantum 54 2.3. Linear optical response of canonical plasmonic nanostructures addressed within TDDFT finite-size effects [ 250 – 252 ] break the invariance with particle size of the spectral profile of σabs ( ω )predicted by classical (nonretarded) LRA calculations (panel b). First, decreasing the size of the MNP within TDDFT produces a redshift of the DP resonance 18 from ωDP ∼ 3 . 25 eV (red, lower spectrum) to ωDP ∼ 3 eV (purple, upper spectrum), in contrast to the size-invariant ωDP = ωp/√3≈ 3 . 25 eV classical LRA value 19 (see Subsection 1.2.3). This redshift is a consequence of the spill out of the induced charges due to nonlocal dynamic screening shown in Figure 2.2b. One can interpret in an intuitive picture that the size of the MNP is effectively increased [ 194 ], which for the same number of electrons reduces the electron density and, thus, also the classical value of ωp [Eq. (1.12) ] and ωDP . The impact of this spill out, and thus the redshift, is more significant the smaller the MNP is [96,256,257]. Further, the width of the DP resonances, κ , obtained with TDDFT increases with decreasing the size of the MNP [ 258 ]. For small-sized MNPs, as considered in Figure 2.3, the width of the plasmon resonance κ obtained within TDDFT is determined by surface-enabled Landau damping [ 64 , 66 , 194 , 259 – 261 ], consisting in the plasmon decay into single electron–hole excitations caused by the scattering at the MNP surface (see also Figure 2.1b). Indeed, the MNP surface provides the momentum required for electrons to excite an intraband transition within the conduction band and create an electron–hole pair [ 262 , 263 ], which is forbidden in the bulk due to momentum conservation. This quantum surface effect gains more importance with decreasing the radius a of the MNP. Indeed, the plamon resonance width associated with Landau damping scales as ∼a−1 [ 264 – 267 ] (see Subsection 2.5.1). Finally, the fragmented shape of σabs ( ω )obtained for the smallest MNP considered in Figure 2.3a (purple, Ne = 338 electrons) is also a consequence of surface-enabled Landau damping. For such small MNPs the energy difference between different single electron–hole transitions is larger than the broadening η = 0 . 07 eV [Eq. (2.56) ] accounting for dissipation processes in the system [ 268 ]. These single-electron features gradually disappear with increasing the number of conduction electrons, because the spectrally close electron–hole transitions merge with each other. In this situation, a Lorentzian resonance profile of σabs is obtained, and thus the system is said to exhibit a “better-developed” plasmonic behavior with increasing size [203,269]. 2.3.2 Dimers of spherical metallic nanoparticles We next consider a dimer composed by two identical spherical metallic nanoparticles (MNPs) of radius a , as schematically shown in Figure 2.4. The electronic structure of the MNPs is described within the jellium model by considering a closed-shell configuration as introduced in Subsection 2.3.1. The gap separation between the 18 This is not the case for noble metals such as Au or Ag, where interband transitions involving d -band electrons give rise to a blueshift of ωDP instead of a redshift with decreasing size of the MNP [80,253–255]. 19 In Figure 2.3b we use ωp = 5 . 63 eV so that the value of ωDP obtained classically coincides with that obtained in Figure 2.3a within TDDFT for the largest MNP. 55 Chapter 2. Quantum many-body description of light–matter interaction with j, m, q the quantum (integer) numbers that define a KS state. The timeindependent KS equation is given in this coordinates by [Eq. (2.5)] ˆ H[n0(ρ)] −1 2 ∂2 ∂z2Ψ0 j,m,q(ρ, z, φ) = ϵj,m +1 22π Lq2!Ψ0 j,m,q(ρ, z, φ), (2.70) where ϵj,m are the eigenvalues of the radial part of the Hamiltonian ˆ H[n0(ρ)], ˆ T+VH[n0(ρ)] + Vxc[n0(ρ)] | {z } ˆ H[n0(ρ)] ψ0 j,m(ρ) = ϵj,mψ0 j,m(ρ),(2.71) and the ρ-space kinetic-energy operator ˆ Tis expressed as ˆ T=−1 21 ρ ∂ ∂ρ ρ∂ ∂ρ−m2 ρ2.(2.72) The equilibrium electron density n0(ρ)is given by n0(ρ)=2 X {j,m,q}∈occ |Ψ0 j,m,q(ρ, z, φ)|2 = 2 X {j,m,q}∈occ 1 2π|ψ0 j,m(ρ)|21 LX q Θ EF− ϵj,m +1 22π Lq2!!, (2.73) where EF is the Fermi energy, Θrepresents the Heaviside step function, and the factor 2 is due to the spin degeneracy. Because of the cylindrical symmetry, the ±mstates are degenerate so that we only solve Eqs. (2.70) and (2.71) for m≥0. The case of the infinite cylindrical nanowire can be obtained by taking L→ ∞ in Eqs. (2.69) - (2.73) . In such a situation, the summation over the quantum number q in Eq. (2.73) can be transformed into an integral over the continuous variable kz=2π Lq, and n0(ρ)can be expressed as n0(ρ) = X {j,m≥0}∈occ χj,m|ψ0 j,m(ρ)|2,(2.74) with the statistical factor χj,m accounting for spin and ±m degeneracy (see Subsection 2.3.1), as well as for the degeneracy due to the electron motion along the z-axis, χj,m =   1 π2p2(EF−ϵj,m)for m= 0 2 π2p2(EF−ϵj,m)for m>0 .(2.75) Similarly to the equilibrium density for spherical MNPs in Subsection 2.3.1 [Eq. (2.48) ], n0 ( ρ )of the cylindrical nanowire, given by Eq. (2.74) , only depends upon the radial coordinate ρ . The orbitals ψ0 j,m ( ρ )of the individual nanowire are therefore obtained from Eqs. (2.71) and (2.74) using a self-consistent procedure 62 2.3. Linear optical response of canonical plasmonic nanostructures addressed within TDDFT based on the KS scheme described in Subsection 2.3.1, where the Hamiltonian ˆ H [ n0 ( ρ )] in Eq. (2.71) is diagonalized by expressing ˆ T [Eq. (2.72) ] in real space with finite differences [75,130]. Once the time-independent KS orbitals ψ0 j,m ( ρ )of the individual nanowire are obtained from Eqs. (2.71) and (2.74) , we use the WPP method described in Subsection 2.2.1 and propagate in real time the KS orbitals ψp subjected to an external time-dependent potential. The time-dependent potential Vext ( x, y, t ) depends on the ( x, y )-coordinates so that the translational invariance of the system along the z -axis is preserved. For the time propagation, we discretize ψp on an equidistant mesh in Cartesian coordinates, ψp≡ψp ( x, y, t ), where the initial conditions ψp ( x, y, t = 0) ≡ψ0 j,m ( ρ, φ )are given by the KS orbitals of the ground state. Note that when expressing the time-dependent KS orbitals in Cartesian coordinates, the quantum number p replaces the pair {j, m} used in cylindrical coordinates for the ground state. One of the advantages of using Cartesian coordinates ψp≡ψp ( x, y, t )for the time propagation is that it is possible to directly apply the same algorithm as for the individual nanowire to study e.g. the optical response of a pair of parallel nanowires (nanowire dimer, Figure 2.6a) under the influence of any external potential Vext ( x, y, t )depending on ( x, y ). This algorithm basically consists in applying the pseudospectral FGH method [ 230 , 236 , 237 ] to calculate the kinetic-energy operator ˆ T as well as the Hartree potential VH[n(x, y, t)] [Eq. (2.23)]. We obtain the optical response of nanowires to a spatially-constant x -polarized external electric field of amplitude E0 , corresponding to the following external potential [see Eq. (2.51)]: Vext(x, y, t) = δ(t)E0∆t x. (2.76) Mesh steps of the order of ∆ x = ∆ y∼ 0 . 5 a0 and a time step of ∆ t∼ 0 . 1 au are typically used in this thesis for the cylindrical geometry. The KS orbitals ψp ( x, y, t ) evolve in time according to the time-dependent KS equations [Eq. (2.19) ], with the time-dependent electron density n(x, y, t)expressed in Cartesian coordinates as n(x, y, t) = X k∈occ χp|ψp(x, y, t)|2.(2.77) In Eq. (2.77), the statistical factors χpare now given by χp=2 πq2(EF−ϵp),(2.78) where ϵpis the eigenenergy of the ground-state KS orbital ψp(x, y, t = 0). Optical response of individual cylindrical metallic nanowires and dimers We show here the TDDFT results for the linear optical response of an individual metallic nanowire and nanowire dimers [ 59 , 119 , 273 – 276 ] consisting of ˜ N = 240 a−1 0 63 Chapter 2. Quantum many-body description of light–matter interaction 0 200 400 600 800 1000 1200 1400 1600 1800 Figure 2.6:(a) Geometry of the individual cylindrical nanowire (left) and the dimer composed by two parallel nanowires (right). The radius of the nanowires is Rc , which is determined by the number of conduction electrons ˜ N per unit length in the z -direction according to Eq. (2.68) . The cylinders are infinite along the z -axis (translationally invariant), and are separated by a gap distance D along the x -axis in the dimer configuration. (b) Absorption cross-section spectrum σabs ( ω ) /L per unit length L in the z -direction of an individual nanowire characterized by a Wigner–Seitz radius rs = 3 . 02 a0 and ˜ N = 240 a−1 0 ( Rc≈ 94 . 3 a0 ). (c) Same as in (b) but for the dimer configuration where the gap separation is varied from D = 14 a0 to D = 40 a0 . An attenuation parameter η= 0.07 eV [Eq. (2.56)] is used to perform the Fourier transform. electrons per unit length in the z -direction (radius Rc≈ 94 . 3 a0 ). We use a jellium model with a Wigner–Seitz radius rs = 3 . 02 a0 ( ωp≈ 8 . 98 eV ) characteristic of the conduction electron density of gold. In addition, we introduce a stabilizing potential 21 [ 216 ] inside the metal so that the work function WF = 5 . 5 eV of gold is retrieved [ 277 ]. Note that, since the jellium model does not account for optical transitions involving localized d -band electrons (see Subsection 2.1.2), the plasmonic response obtained here presents important differences with respect to the results that would be obtained for a more exact model of gold. However, this description of the metallic nanowires still allows one to predict electron transport properties of gold junctions in off-resonant excitation conditions, as implemented in recent works [ 77 , 210 , 278 ] when studying the electron-currents dynamics induced by ultrafast 21 Applying a stabilizing potential simply consists in introducing a constant potential inside the metal in the effective potential given by Eq. (2.17) both in the ground-state and in the time-dependent calculations. This procedure does not introduce any additional computational difficulty. 64 2.4. Nonlinear effects in the optical response of spherical plasmonic nanoparticles addressed within TDDFT electromagnetic fields in plasmonic gaps. Here we are interested in the general trends of the optical response of these systems, where similar physical effects as described in Subsections 2.3.1 and 2.3.2 for spherical MNPs are also expected to be present. Figure 2.6b shows the absorption cross-section spectrum σabs ( ω ) /L per unit length L in the z -direction of an individual metallic nanowire of radius Rc≈ 94 . 3 a0 ( ≈ 5 nm ). Due to the relatively large size of the nanowire, σabs ( ω ) /L does not exhibit single electron–hole transition features and the plasmonic response is well developed. A single peak associated with the DP resonance ( m = 1) of the nanowire emerges at ωDP = 6 . 24 eV , slightly below the classical nonretarded prediction, ωDP = ωSP = ωp/√2 = 6 . 35 eV22 [ 22 ], as a consequence of the spill out of the induced charges (see Subsection 2.3.1). On the other hand, Figure 2.6c displays the absorption spectrum σabs ( ω ) /L of a dimer formed by two parallel nanowires, identical to that discussed in Figure 2.6b. The gap separation, D , ranges from D = 14 a0 to D = 40 a0 ( D∼ 0 . 75 − 2 . 1 nm ). The overall qualitative behavior of σabs ( ω ) /L with reducing gap distance D is almost identical to either the classical LRA (Figure 1.5a) or TDDFT (Figure 2.5c) results of the spherical MNP dimers, and can be understood using the electromagnetic hybridization picture introduced in Subsection 1.2.3 [ 150 , 151 ]. Note that, due to the azimuthal symmetry of the individual nanowires, here the magnetic quantum number m plays a similar role as the multipole order ℓ in spherical MNPs. For the largest gap separation, D = 40 a0 ( D∼ 2 . 1 nm ), two distinct modes emerges: a BDP at ωDP ∼ 5 . 35 eV formed from the hybridization of the DP modes ( m = 1) of the individual nanowires, and a broad pseudomode at ωPSM ∼ 6 . 2 eV formed by the hybridization of nearly-degenerated higher-order modes ( m > 1). This pseudomode is slightly redshifted with respect to the surface plasmon frequency. For smaller gap separations, D∼ 28 a0 , σabs ( ω ) /L exhibits three well-defined peaks because another distinct resonance corresponding to the BQP at ωBQP ∼ 5 . 9 eV emerges. The BDP and the BQP resonances redshift with reducing D because of the increasing attractive interaction between the charges of opposite sign across the junction [126]. 2.4 Nonlinear effects in the optical response of spherical plasmonic nanoparticles addressed within TDDFT In previous sections, we focused on the linear optical response of MNPs, which are shown to strongly enhance the strength of the incident electric field at optical frequencies because of the excitation of plasmon resonances. However, plasmonic nanostructures also exhibit an efficient nonlinear optical response to 22 While the DP resonance in spherical MNPs ( ℓ = 1) is classically at ωDP = ωℓ=1 = ωp/√3 within the nonretarded approximation, the DP resonance in cylindrical nanowires ( m = 1) is at the surface plasmon frequency, ωDP =ωm=1 =ωSP =ωp/√2[22]. 65 Chapter 2. Quantum many-body description of light–matter interaction strong illumination intensities, thus being good candidates for the fabrication of nanodevices based on nonlinear optics [ 36 ]. One manner to theoretically study the nonlinear optical response of MNPs is to adopt the classical scheme introduced in Section 1.3 with the use of the nonlinear hyperpolarizabilities α(n) ,ad hoc parameters that characterize the nonlinear optical response of the system. In this section, we show that the real-time TDDFT approach based on the WPP method employed in previous sections to study the linear optical response of plasmonic structures can also be used to analyze nonlinear effects in the optical response without using ad hoc parameters. As an example, we consider the nonlinear optical response of an individual spherical MNP. We describe the electronic structure of the MNP within the jellium model introduced in Subsection 2.1.2 using a Wigner–Seitz radius rs = 4 a0 of sodium. We address a spherical MNP that contains 1074 conduction electrons, with a radius a = 40 . 96 a0 ( ≈ 2 . 2nm). The numerical implementation is almost identical to the one used in Subsection 2.3.1. The only difference is that for this study of nonlinear effects, instead of using the external potential Vext ( r, t )given by Eq. (2.51), here we use the following one: Vext(r, t) = E0rcos θcos(ω(t−t0)) e−(t−t0 σ)2 ,(2.79) which corresponds to the potential experienced by an electron interacting with an incident Gaussian laser pulse polarized along the z -axis. The fundamental frequency of the external illumination ω = 1 . 585 eV in this section is half of the DP frequency ωDP = 3 . 17 eV (see Figure 2.3a), the duration of the pulse is σ = 5 × 2 π/ω , and the arrival time of the pulse t0is t0= 5σ. In order to analyze the nonlinear optical response of the individual MNP, we calculate within TDDFT the time evolution of the induced electron density δn ( r, t ) and obtain the time-dependent induced dipole moment, p ( t )[Eq. (2.33) ], as well as the electric near field (induced field), Eind ( r, t )[Eq. (2.34) ], created by the MNP in response to an incident electromagnetic pulse. The frequency-resolved quantities are then obtained from the time-to-frequency Fourier transform given by Eq. (2.36) , δn(r,Ω) = Zdt δn(r, t)eiΩte−(t−t0 σ)2 , p(Ω) = Zdt p(t)eiΩte−(t−t0 σ)2 , Eind(r,Ω) = Zdt Eind(r, t)eiΩte−(t−t0 σ)2 , (2.80) where in this case the filter F(t)is given by F(t) = e−(t−t0 σ)2 .(2.81) The Gaussian filter F ( t )introduced in Eq. (2.81) partially accounts for decay and dephasing processes of the collective density oscillations that are not included in the present ALDA-TDDFT approach [ 72 , 241 , 242 ] (see Section 2.2), and allows us 66 2.4. Nonlinear effects in the optical response of spherical plasmonic nanoparticles addressed within TDDFT Figure 2.7: Nonlinear optical response of the individual spherical MNP as calculated within TDDFT for an incident z -polarized Gaussian electromagnetic pulse with fundamental frequency ω = 1 . 585 eV (half of the frequency of the MNP dipolar plasmon ωDP = 3 . 17 eV ).The intensity is I0 = 10 8 W cm −2 (dashed blue line) or I0 = 10 10 W cm −2 (red line). Panel (a) shows the square of the induced dipole moment |p (Ω) |2 , and panel (b) the absolute value of the spectrum of the electric near field |Eind ( r, Ω) | induced at the z -axis at 18 a0 ( ≈ 0 . 95 nm ) from the MNP surface. The symbol Ωrepresents the frequency of the induced electromagnetic fields in response to the incident illumination with fundamental frequency ω. to reach convergent spectral response at high-harmonic frequencies. This approach is justified because the fundamental frequency is strongly detuned from the DP resonance of the MNP, so that no electron-density oscillation and high-harmonic generation is expected when the laser is switched off. Consistently, we apply a Gaussian filter given by the envelope of the incident pulse. Moreover, in Eq. (2.80) we use the symbol Ω(compare Eq. (2.80) with Eq. (2.36) ) to refer to the oscillation frequency of the electromagnetic fields (and dipole moments) induced by the MNP due to the nonlinear optical response to the external excitation oscillating at ω. The nonlinear optical response of the individual spherical MNP is displayed in Figure 2.7. In panel (a), we show the intensity spectrum of the induced dipole moment |p (Ω) |2 , which is proportional to the power of light emitted to the far field [Eq. (1.46) ]. In panel (b), we show the spectrum of the induced near field |Eind ( r, Ω) | at the z -axis, at 18 a0 ( ≈ 0 . 95 nm ) from the MNP surface. Results are obtained for an incident Gaussian electromagnetic pulse with intensity I0 = 10 8 W cm −2 ( E0 = 4 . 8 × 10 −5au , dashed blue line) and I0 = 10 10 W cm −2 ( E0 = 4 . 8 × 10 −4au , red line), averaged over the duration of the pulse σ . The corresponding energy per incident pulse is well below the documented damage threshold of small MNPs [ 279 – 281 ]. The induced dipole moment |p (Ω) |2 in Figure 2.7a exhibits only odd harmonics n = 1 , 3 , 5 . . . . Thus, only odd multiples Ω = ω, 3 ω, 5 ω, . . . of the incoming frequency ω are emitted by the system into the far field, consistently with the inversion symmetry of the MNP that prevents even-harmonic generation [ 152 , 282 ] (see Section 1.3). Overall, the nonlinear response |p (Ω) |2 for I0 = 10 10 W cm −2 is several orders of magnitude larger than that for I0 = 10 8 W cm −2 . This large increase is in accordance with the In 0 dependence of |p (Ω = nω ) |2 , expected from the standard theory of nonlinear optics as described in Section 1.3 [152]. In contrast to the far-field response, both odd and even harmonics are present in 67 Chapter 2. Quantum many-body description of light–matter interaction Figure 2.8: Color maps of the real part of the induced electron density δn ( r, Ω) (left), of the radial component of the electric near field Eind ( r, Ω) (center), and of the tangential component of Eind ( r, Ω) (right) induced at the fundamental, second, third, and fourth-harmonic frequency by a z -polarized Gaussian electromagnetic pulse with fundamental frequency ω = 1 . 585 eV and intensity I0 = 10 10 W cm −2 incident at the individual spherical MNP. Results are rotationally symmetric with respect to the z -axis, and they are shown in the ( x, z )-plane normalized to unity. the spectrum of the electric near field induced by the individual MNP (Figure 2.7b). Indeed, at the metal–vacuum interface the inversion symmetry is locally broken, and short-range even-harmonic electric fields can be induced close to the MNP surface [ 283 – 287 ]. As expected, the |Eind ( r, Ω) | is orders of magnitude larger for I0= 1010 W cm−2(red line) than for I0= 108W cm−2(blue). The color maps of the induced electron density δn ( r, Ω) and of the electric near field Eind ( r, Ω) induced by the incident z -polarized Gaussian electromagnetic pulse are shown in Figure 2.8 for the fundamental, second, third, and fourth harmonics. The induced charge density δn ( r, Ω = nω )of the n -th harmonic and the corresponding near field Eind ( r, Ω = nω )are shown in the ( x, z )-plane. Because of the symmetry of the configuration, the calculated color maps are 68 2.5. Semiclassical surface-response formalism (SRF) independent of a rotation around the z -axis. At odd harmonics ( n = 1 , 3), the induced charge densities are antisymmetric with respect to the ( x, y )-plane, δn ( x, y, z, nω ) = −δn ( x, y, −z, nω ), which results in a net dipole moment (see Figure 2.7a). In contrast, a quadrupolar-like near field and symmetric chargedensity, δn ( r, nω ) = δn ( -r, nω ), are induced at even harmonics ( n = 2 , 4) [ 288 – 294 ]. The dipole moment p (Ω) is zero in this case, and thus there is no emission into the far field at even harmonics. Thus, despite the second harmonic being at resonance with the dipolar plasmon of the MNP, the latter can not be excited because of the symmetry selection rules. Another consequence of the symmetry selection rules is that, for any point located in the ( x, y )-plane, the even-harmonic near field in that point is oriented perpendicularly to the z -polarized incident pulse. We use in Chapter 6the insights obtained in this section to study nonlinear effects when the spherical MNP is coupled to a QE located nearby, demonstrating that the presence of the QE enables the emission into the far field at the second-harmonic frequency, otherwise forbidden because of the inversion symmetry of the MNP. 2.5 Semiclassical surface-response formalism (SRF) In this section, we introduce the semiclassical surface-response formalism (SRF) employed in this thesis, which allows one to account for quantum surface effects using much less computationally-demanding calculations as compared to the TDDFT methodology. In brief, the SRF is an extension of the classical LRA (Chapter 1) that incorporates surface-response corrections at the metal–dielectric interface in the solution of Maxwell’s equations (Section 1.1) by means of the so-called Feibelman parameters obtained from quantum-mechanical calculations [ 95 – 98 ]. These parameters, commonly denoted as d⊥ and d∥ , were first introduced by Peter Feibelman in the 1980s [ 98 ], and have received renewed attention during the last years [ 97 , 124 , 129 , 276 , 295 – 302 ] due to their usefulness to study systems that have recently become experimentally feasible. Indeed, the semiclassical SRF based on the Feibelman parameters allows us to account for the gradual variation of the induced electron density across the metal–vacuum interface (see Figure 2.2b), in contrast to the classical LRA described in Chapter 1that considers the polarization charges to be located strictly at the metal boundary of infinitesimal width. The Feibelman parameters d⊥ and d∥ are usually defined by considering a semi-infinite metal surface [ 54 , 55 , 62 , 98 ], although the expressions for d⊥ and d∥ have been also proposed for other geometries such as spherical MNPs [ 95 ] 23 . To our knowledge, the Feibelman parameters have been so far computed in the literature within the long-wavelength approximation, which consists in neglecting the nonlocality of the optical response in the direction parallel to the metal– dielectric interface. Thus, within the long-wavelength approximation, d⊥≡d⊥ ( ω ) and d∥≡d∥ ( ω )solely depend on the excitation frequency, ω , and not on the 23 In Chapter 4we analyze the Feibelman parameter d⊥for a cylindrical metallic nanowire. 69 Chapter 2. Quantum many-body description of light–matter interaction wavenumber parallel to the metal surface, k∥ . In this thesis we thus refer to the k∥-independent d⊥(ω)and d∥(ω)as the nondispersive Feibelman parameters. We adopt the long-wavelength approximation in this section and in Chapter 3, where we identify situations where this approximation becomes inaccurate. In Chapter 4, we calculate the Feibelman parameters as a function of both ω and k∥ and propose a dispersive SRF that overcomes the shortcomings of the long-wavelength approximation. The nondispersive Feibelman paramaters defined for a semi-infinite planar metal– vacuum interface at z = 0 (Figure 1.3a) are usually expressed as [ 62 , 81 , 124 , 296 ]: d⊥(ω) = Rdz z δn(z, ω) Rdz δn(z, ω),(2.82a) d∥(ω) = Rdz z ∂ ∂z j∥(z, ω) Rdz ∂ ∂z j∥(z, ω),(2.82b) with δn ( z, ω )the quantum-mechanical (complex-valued) electron density induced in response to the external excitation, and j∥ ( z, ω )the parallel-to-the-surface component of the associated induced electron current density. Other definitions equivalent to Eq. (2.82) that express the parameters in terms of the electromagnetic fields and nonlocal dielectric functions have been also used [ 62 , 98 , 296 , 303 ]. In Eq. (2.82) , the real part of d⊥ ( ω )( Re{d⊥ ( ω ) } ) corresponds to the position of the centroid of the induced charge density with respect to the positive background edge of the metals (see schematic representation in Figure 2.9a), while the imaginary part ( Im{d⊥ ( ω ) } ) is related to surface-enabled Landau damping [ 62 ] (see below). On the other hand, Re{d∥ ( ω ) } represents the position of the centroid of the normal derivative of the electron current parallel to the metal surface. The Feibelman parameters d⊥ ( ω )and d∥ ( ω )are surface-response functions inherent to a specific metal, but also dependent on the surrounding material [ 62 ]. Importantly, d∥ ( ω )given by Eq. (2.82b) vanishes for charge-neutral planar surfaces [ 62 , 304 ], and it is also expected to be much less important than d⊥ ( ω ) for curved surfaces [ 95 ]. In this thesis we thus consider d∥ ( ω ) = 0, consistent with the approximation adopted in other studies [ 81 , 124 , 296 , 300 ]. We show in Figure 2.9b the nondispersive Feibelman parameter d⊥ ( ω )used in this section and in Chapter 3, obtained in ref. 296 for a sodium planar surface ( rs = 4 a0 ) surrounded by vacuum 24 . The Feibelman parameter d⊥ ( ω )in Figure 2.9b shows a resonance at ω∼ 4 . 7 eV , associated with the excitation of the Bennet plasmon at ω∼ 0 . 8 ωp [ 79 , 305 ], also referred to as the multipole surface plasmon [ 62 , 306 ] (not to be confused with localized multipole plasmons supported e.g. by spherical MNPs). Kramers-Kronig relations connect the real (blue line) and imaginary (red line) parts of d⊥(ω)[307]. The advantage of the SRF is that, once the Feibelman parameters are obtained 24 Ref. 124 provide the parametrization to express the data of d⊥ ( ω )obtained in ref. 296 within TDDFT as a sum of Lorentzian functions. 70 2.5. Semiclassical surface-response formalism (SRF) -6 -4 -2 0 2 4 6 8 10 Figure 2.9:(a) Schematic representation of the nondispersive Feibelman parameter d⊥ ( ω )in a planar metal–vacuum interface. Re{d⊥ ( ω ) } corresponds to the position of the centroid of the induced charge δn ( z, ω )with respect to the classical metal surface located at z = 0.(b) Real (blue) and imaginary (red) parts of the nondispersive Feibelman parameter d⊥ ( ω )obtained in ref. 296 for a planar sodium–vacuum interface. The parameter shown in panel (b) is used in this section and in Chapter 3, where we study quantum surface phenomena in the interaction between QEs and MNPs. for a planar surface of a given material using quantum-mechanical methods, these parameters can in principle be applied to account semiclassically for quantum surface effects in the optical response of arbitrary-shaped MNPs. Within the semiclassical SRF, the electromagnetic problem is addressed by solving Maxwell’s equations [Eq. (1.5) ] with the use of local dielectric functions (e.g. a Drude dielectric function, Eq. (1.10) ), and introducing a set of modified boundary conditions at the metal–dielectric interfaces that differ from those used within the classical LRA [Eqs. (1.16) and (1.18) ]. d⊥ ( ω )and d∥ ( ω )can be related, respectively, to a surface polarization oriented perpendicularly to the interface and to a parallel surface current [ 97 , 296 ] leading to the following modified boundary conditions (see also ref. 124): ˆ n×ESRF out −ESRF in =−d⊥(ω)ˆ n×∇hˆ n·ESRF out −ESRF in i,(2.83a) ˆ n·DSRF out −DSRF in =d∥(ω)∇·hˆ n×DSRF out −DSRF in ׈ ni,(2.83b) ˆ n×HSRF out −HSRF in =iωd∥(ω)hˆ n×DSRF out −DSRF in ׈ ni,(2.83c) ˆ n·BSRF out −BSRF in = 0,(2.83d) where ˆ n is the normal unit vector pointing outwards from the metal boundary, and the superscript “SRF” denotes that the fields are calculated within the semiclassical SRF. Note that, as mentioned above, we consider d∥ ( ω )=0, so that in practice we only consider the modifications introduced by Eq. (2.83a) into the boundary conditions. The rest of the boundary conditions remain the same as in the classical 71 Chapter 3. Quantum surface effects in the electromagnetic coupling between quantum emitters and metallic nanoparticles which in the weak-coupling regime determines the Purcell factor [ 101 , 102 ] and Lamb shift [ 103 , 104 ] in the emission by the QE produced by the plasmonic environment. Typically, a local dielectric function obtained experimentally or from simple theoretical approaches, such as the Drude model (Subsection 1.1.1), can be used to characterize the MNPs. However, such a dielectric function does not account for quantum phenomena relevant in MNPs of small characteristic dimensions such as electron spill-out, surface-enabled Landau damping, or nonlocal dynamical screening [ 10 , 61 , 64 , 66 , 96 , 126 , 209 , 212 , 219 ], introduced in Section 2.3 when analyzing the linear optical response of MNPs. These nonclassical phenomena, inherent to the quantum nature of electrons in metals, are also expected to influence the QE–MNPs electromagnetic interaction when small MNPs and QE–MNPs distances are considered [65,124,300,330–336]. In this chapter, we use the TDDFT approach [ 58 , 79 , 80 , 202 , 210 , 214 , 337 – 341 ] based on the WPP method introduced in Chapter 2to provide a fundamental description of the electromagnetic coupling between QEs and canonical MNPs. The QE–MNP separations are set sufficiently large so that electron tunneling is negligible and thus the electromagnetic interaction mainly determines the optical properties of the system. However, we consider QE–MNP separations small enough for the aforementioned quantum surface effects to be important. We first use TDDFT to calculate the self-interaction Green’s function ˆ G ( rQE,rQE, ω )that governs the QE–MNPs coupling. A comparison with classical LRA results (Subsection 1.1.1) reveals the importance of quantum effects. Further, we compare TDDFT results and the semiclassical nondispersive SRF (Section 2.5) that incorporates quantum surface-response corrections via the Feibelman parameter d⊥ ( ω )obtained in the long-wavelength approximation, which allows us to identify surface-enabled Landau damping and spill-out of the induced electron density as the dominant quantum mechanisms dramatically influencing the electromagnetic QE–MNPs interaction. TDDFT also provides a benchmark to establish the validity range of the (standard) nondispersive SRF to adequately account for the dominant quantum phenomena arising in the electromagnetic interaction between QEs and MNPs. The nondispersive SRF, as used in this chapter, neglects the nonlocal optical response in the direction parallel to the metal surface (long-wavelength limit), enabling a very efficient implementation of nonlocality in nanoscale geometries [ 97 ]. In this chapter, we identify situations of very small QE–MNP distances where the nondispersive implementation of the SRF eventually fails, indicating that the dispersion of the Feibelman parameters with respect to the wavenumber parallel to the surface needs to be considered (see Chapter 4). 3.1 System and methods We analyze the electromagnetic coupling between a point-like QE and two different canonical plasmonic nanostructures. First, we consider in Subsections 3.2.1,3.2.2, and 3.2.3 the case of an individual spherical MNP, where the QE is placed at a distance d from the MNP surface (Figure 3.1a). Then, in Subsection 3.2.4, we 78 3.1. System and methods Figure 3.1: Sketch of the systems studied in this chapter, consisting of a QE modeled as a point dipole pd placed (a) at a distance d from the surface of an individual spherical MNP, and (b) at the center of a gap of size D formed by two identical spherical MNPs. The point dipole is oriented along the z -axis, which is also the axis of the dimer. Each MNP is represented within the free-electron jellium model using a Wigner–Seitz radius rs = 4 a0 that corresponds to sodium, and contains Ne= 4458 conduction electrons resulting in a radius a= 65.83 a0(≈3.5nm). analyze a dimer composed by two identical spherical MNPs separated by a gap distance D (Figure 3.1b), with the QE situated at the center of the gap ( z = 0). We define the coordinates in the same way as in Subsections 2.3.1 and 2.3.2, such that the center of the nanoparticle(s) and the QE is at the z -axis. The entire system is surrounded by vacuum. We consider closed-shell jellium MNPs of radius a = 65 . 83 a0 ( ≈ 3 . 5 nm ), which results in a well-developed plasmonic response [ 203 ] in the TDDFT simulations. The surface-to-emitter distance is sufficiently large to ensure that the electron densities of the MNPs at the position of the QE are negligible, and therefore there is no electron tunneling [ 43 , 67 , 68 , 83 , 86 , 342 ]. Specifically, we consider surface-to-emitter distances d in the range of d = 10 − 42 a0 ( ≈ 0 . 5 − 2 . 2 nm ) for the case of the individual MNP, and gap separations of D = 2 d = 20 − 45 a0 ( ≈ 1 . 1 − 2 . 4 nm ) for the dimer structure. We use the point-dipole approximation described in Subsection 1.4.1 to model the QE and thus neglect its spatial extent. We focus on the study of the self-interaction Green’s function ˆ G ( rQE,rQE, ω ) defined according to Eq. (1.56) , which provides the electric field Eself ( rQE, ω ) created by the metallic nanostructure at a position rQE in response to a point dipole pd located at the same position and oscillating at a frequency ω . Importantly, ˆ G ( rQE,rQE, ω )determines the total decay rate (Γ) and the Lamb shift (∆ ωQE ) of a QE that interacts weakly with the plasmonic nanostructure [104,170–172], Γ = γ0+γnr QE + 2|µQE|2Im{ˆ k·ˆ G(rQE,rQE, ω =ωQE)·ˆ k},(3.1a) ∆ωQE =−|µQE|2Re{ˆ k·ˆ G(rQE,rQE, ω =ωQE)·ˆ k},(3.1b) with γ0 , γnr QE and µQE the spontaneous decay rate in vacuum, the non-radiative 79 Chapter 3. Quantum surface effects in the electromagnetic coupling between quantum emitters and metallic nanoparticles intrinsic loss rate, and the transition dipole moment of the QE along the ˆ k -direction, respectively [ 22 ]. The self-interaction Green’s function ˆ G ( rQE,rQE, ω )in Eq. (3.1a) and Eq. (3.1b) is evaluated at the transition frequency of the QE, ω = ωQE . The unit vector ˆ k defines the orientation of pd . We note that the enhancement of the decay rate due to the QE–MNPs coupling is often normalized by γ0 , which gives the Purcell factor FP=Γ−γnr QE γ0 = 1 + 3c3 2ω3 QE Im{ˆ k·ˆ G(rQE,rQE, ω =ωQE)·ˆ k},(3.2) cbeing the speed of light in vacuum [102,343]. We restrict our analysis to the case of a point dipole oriented along the z -axis, pd = pdˆ z , where ˆ z is the unit vector along the z -axis. Moreover, we use rQE = rQE ˆ z , so that ˆ G ( rQE,rQE, ω )can be considered as a scalar, ˆ G ( rQE,rQE, ω ) ≡G ( rQE,rQE, ω ), because of the symmetry of our system. In the upcoming subsections, we briefly summarize the key aspects of the three different models used in this chapter to obtain G ( rQE,rQE, ω ), namely the TDDFT, the classical LRA, and the semiclassical SRF based on the nondispersive Feibelman parameters. Further details on the methodologies adopted in this chapter are provided in Chapters 1and 2. 3.1.1 Time-dependent density functional theory (TDDFT) In the TDDFT calculations we describe the MNPs within the jellium model of free-electron metals [ 194 , 200 ] as introduced in Subsection 2.1.2. We use a Wigner– Seitz radius equal to that of sodium, rs = 4 a0 , which allows us to use the values of the nondispersive Feibelman parameters obtained in ref. 296 (Figure 2.9b) to calculate G ( rQE,rQE, ω )within the SRF and make the comparison between the TDDFT and SRF results. Note also that the DP resonance of sodium MNPs lies at optical frequencies ( ωDP ∼ 3 eV , see Figure 2.3), very close to that of gold MNPs, thus placing the results within the frequency range relevant for actual applications in plasmonics. The closed-shell MNPs contain Ne = 4458 conduction electrons, which sets the radius of the background jellium edge to a= 65.83 a0(≈3.5nm). In order to obtain the self-interaction Green’s function within TDDFT, GTDDFT ( rQE,rQE, ω ), we use the WPP method explained in Subsection 2.2.1 to solve the time-dependent KS equations given by Eqs. (2.19) and (2.20) in time domain, and to calculate the time evolution of the electron density n ( r, t ). Instead of Eq. (2.51) used to model plane-wave illumination in Chapter 2, the external potential Vext ( r, t )that drives the system [Eq. (2.21) ] is given in this chapter by an impulsive potential created by a point dipole, Vext(r, t) = −pd∆tˆ z·r−rQE |r−rQE|3δ(t),(3.3) with δ ( t )the Dirac delta function and the amplitude pd sufficiently small to ensure a linear response. Vext ( r, t )appears with a minus sign in Eq. (3.3) because it is 80 3.1. System and methods acting on electrons with charge qe = − 1. For the individual MNP, we employ the method described in Subsection 2.3.1 and represent the occupied KS orbitals Ψ j ( r, t )using a meshgrid in spherical coordinates, whereas for the dimer we adopt cylindrical coordinates as introduced in Subsection 2.3.2. The propagation time-step is ∆t∼0.05 au. The Hartree potential VH [ n ( r, t )] [Eq. (2.23) ] is calculated for the individual spherical MNP by expressing VH [ n ( r, t )] as a sum of Legendre polynomials Pℓ ( cosθ ), VH[n(r, t)] = ∞ X ℓ=0 1 rVℓ(r, t)Pℓ(cosθ),(3.4) since the system subjected to the external potential Vext ( r, t )given by Eq. (3.3) possesses rotational symmetry with respect to the z -axis, and therefore the response is independent of the azimuthal angle φ . Using the Laplace operator ∇2 in spherical coordinates [Eq. (1.24) ], Poisson’s equation given by Eq. (2.23) can be expressed as: ∞ X ℓ=0 1 rd2 dr2−ℓ(ℓ+ 1) r2 | {z } Aℓ Vℓ(r, t)Pℓ(cosθ) = −4π∞ X ℓ=0 nℓ(r, t)Pℓ(cosθ),(3.5) where nℓ(r, t) = 2 2ℓ+ 1 Z1 −1 (n(r, t)−n+)Pℓ(cosθ)d(cos θ),(3.6) and we express the operator Aℓ = d2 dr2−ℓ(ℓ+1) r2 in matrix form using the method of the Fourier grid Hamiltonian with sine basis functions [ 344 ]. From Eq. (3.5) , Vℓ(r, t)can then be directly obtained from Vℓ(r, t) = −4πA−1 ℓr nℓ(r, t).(3.7) For the dimer geometry we obtain VH [ n ( r, t )] following the procedure described in Subsection 2.3.2 [75]. Finally, we calculate the time-dependent electric field Eind ( r = rQE, t )induced by the metallic nanostructure at the position rQE of the QE using VH [ n ( r, t )] according to Eq. (2.34) . The time-to-frequency Fourier transform [Eq. (2.36) ] finally leads to the frequency-resolved self-interaction Green’s function [compare to Eq. (1.56)], GTDDFT(rQE,rQE, ω) = 1 pd∆tZTf 0 dt Eind(rQE, t)ei(ω+iη/2)t | {z } Eself(rQE,ω) ,(3.8) where Tf = 3500 au is the total propagation time used in our simulations (enough to achieve convergence), and η = 0 . 07 eV accounts for relaxation processes beyond 81 Chapter 3. Quantum surface effects in the electromagnetic coupling between quantum emitters and metallic nanoparticles the ALDA-TDDFT description [ 241 , 242 ] of the many-body dynamics such as the interaction of excited electrons with phonons and many-body inelastic electron– electron scattering events (see Sections 2.2 and 2.3). 3.1.2 Classical local-response approximation (LRA) The classical (nonretarded) local self-interaction Green’s function, GLRA ( rQE,rQE, ω ), is obtained in the frequency domain using the methodology described in Subsection 1.1.3 for the external potential given by Eq. (1.53) , ϕQE ( r ) = pQEˆ z·r−rQE |r−rQE|3 . The specific details to obtain GLRA ( rQE,rQE, ω )for both the individual spherical MNP and the dimer are given in Subsection 1.4.2. In this chapter, we use a Drude-type local dielectric function ε ( ω )to characterize the metal [Eq. (1.10) ], with a plasma frequency ωp = q3 r3 s = 5 . 89 eV and intrinsic damping parameter γp = 0 . 1 eV . The value of parameter γp is obtained from the comparison of the absorption spectrum σabs ( ω )of the individual spherical MNP calculated with TDDFT and with the SRF, as detailed in the following Subsection 3.1.3. In the case of a point-like QE at position rQE oriented in the radial direction and exciting the individual MNP of radius a , the self-interaction Green’s function is given by Eq. (1.58a) [145] GLRA(rQE,rQE, ω) = ∞ X ℓ=1 (ℓ+ 1)2a2ℓ+1 R2ℓ+4 ε(ω)−1 ε(ω) + ℓ+1 ℓ ,(3.9) with R the distance between the position rQE of the QE and the center of the MNP, and ℓthe multipole order of the plasmonic resonance. 3.1.3 Semiclassical surface-response formalism (SRF) The semiclassical SRF employed in this chapter is based on the nondispersive Feibelman parameter d⊥ ( ω )given by Eq. (2.82a) [ 95 , 98 ], a (frequency-dependent) complex-valued function that allows for incorporating quantum surface effects into an otherwise classical description. Indeed, Re{d⊥ ( ω ) } determines the position of the centroid of the induced charge density with respect to the positive jellium edge of the metal, and Im{d⊥ ( ω ) } accounts for surface-enabled Landau damping [ 62 ]. As explained in Section 2.5, in this chapter we consider d∥ ( ω ) = 0 [Eq. (2.82b) ] and use the nondispersive parameter d⊥ ( ω )obtained by Christensen et al. [ 296 ] within the jellium model (Figure 2.9b) for a semi-infinite planar metal surface with rs = 4 a0 (see also ref. 124 ). Using the parameter d⊥ ( ω )obtained for planar surfaces is a reasonable approximation when the radius of curvature of the nanostructure or the typical length of the variation of the external potential along the surface is much larger than the Feibelman parameter. One of the objectives of this chapter is indeed to test the validity of this approximation for situations where the radius of curvature of the system is small and the external potential along the metal surface 82 3.2. Results and discussion (here created by the point-dipole QE) varies rapidly. Similarly to the classical LRA, within the semiclassical SRF one can obtain an analytical solution of the nonretarded self-interaction Green’s function in the present system, GSRF ( rQE,rQE, ω ). The expression for GSRF ( rQE,rQE, ω )can be obtained from Eq. (2.86) and Eq. (2.87) considering the external potential ϕext ( r, ω ) ≡ϕQE ( r )of a point dipole given by Eq. (1.53) , which results in (for d∥(ω)=0) [95,124]: GSRF(rQE,rQE, ω) = ∞ X ℓ=1 (ℓ+ 1)2a2ℓ+1 R2ℓ+4 (ε(ω)−1)1 + ℓ ad⊥(ω) ε(ω) + ℓ+1 ℓ−(ε(ω)−1)ℓ+1 ad⊥(ω).(3.10) For d⊥ ( ω ) /a → 0, i.e., for situations where the radius of the MNP a is much larger than the surface-response correction d⊥ ( ω ),Eq. (3.10) reduces to the classical LRA expression of GLRA(rQE,rQE, ω)given by Eq. (3.9). The value of the damping parameter γp used in the Drude-type dielectric function ϵ ( ω )[Eq. (1.10) ] for the classical LRA as well as for the SRF results is obtained from the fitting of the absorption cross-section spectrum σabs ( ω )of the individual MNP calculated within TDDFT to the result obtained within the nondispersive SRF. The absorption cross section σabs ( ω )within TDDFT is calculated using the methodology described in Subsection 2.3.1 for plane-wave illumination (given by Eq. (2.51) ), and the SRF value is obtained using the polarizability αSRF ( ω )given by Eq. (2.88) . Figure 3.2 shows very good agreement between the TDDFT results (solid blue line) and the nondispersive SRF (dashed red line) when using γp = 0 . 1 eV , thus justifying the value of γp used in this chapter. For simplicity, in the following we generically use G ( rQE,rQE, ω )to refer to any of GTDDFT ( rQE,rQE, ω ), GLRA ( rQE,rQE, ω ), and GSRF ( rQE,rQE, ω ). The methodology used to calculate G ( rQE,rQE, ω )will be clear in the context of each subsection. 3.2 Results and discussion In this section, we present the results for the electromagnetic coupling between a QE and spherical metallic nanostructures. First, in Subsection 3.2.1, we focus on the role of quantum phenomena by comparing quantum TDDFT and classical LRA results for a QE in proximity to an individual spherical MNP. Then, in Subsection 3.2.2, we compare results of TDDFT and the semiclassical nondispersive SRF for the same system. In Subsection 3.2.3, we use the nondispersive SRF to analyze the origin of the observed quantum effects, and discuss the validity range and shortcomings of the long-wavelength Feibelman paramater d⊥ ( ω )that neglects the nonlocality of the optical response in the direction parallel to the metal surface. Finally, in Subsection 3.2.4, we extend the analysis to the case of a dimer of two identical spherical MNPs, showing that the quantum surface effects observed for the individual MNP are also manifested in the dimer configuration and that the 83 Chapter 3. Quantum surface effects in the electromagnetic coupling between quantum emitters and metallic nanoparticles 0 20 40 60 80 100 120 140 160 180 1 1.5 2 2.5 3 3.5 4 4.5 5 TDDFT Figure 3.2: Comparison between the absorption cross-section spectra σabs obtained from TDDFT calculations within the jellium model using a Wigner–Seitz radius rs = 4 a0 (solid blue line), and from the nondispersive SRF using a plasma frequency ωp = 5 . 89 eV and intrinsic damping parameter γp = 0 . 1 eV in the Drude dielectric function given by Eq. (1.10) (dashed red line). The nondispersive SRF results are obtained from Eq. (2.88) , and the TDDFT results are obtained following the procedure described in Subsection 2.3.1. An attenuation factor η = 0 . 07 eV [Eq. (2.56)] is used to perform the Fourier transform within TDDFT. nondispersive SRF does not describe accurately situations of gap distances narrower than D∼1.5nm for the present dimer geometry. 3.2.1 Quantum TDDFT vs. classical LRA We first analyze the quantum effects that influence the self-interaction Green’s function G ( rQE,rQE, ω )obtained for a QE placed in front of an individual spherical MNP. The QE is oriented in the radial direction perpendicular to the MNP surface (see sketch in Figure 3.1a). To identify the quantum effects, we first compare in Figure 3.3 the classical LRA (panels a,b) and the TDDFT (panels c,d) results. We plot both the imaginary (a,c) and real (b,d) parts of G ( rQE,rQE, ω )which are related to the Purcell factor and Lamb shift, respectively [Eqs. (3.1a) and (3.1b) ]. Results are shown as a function of the oscillation frequency of the QE, ω , and the distance dbetween the QE and the surface of the spherical MNP, d=R−a. The classical LRA calculations predict a dependence of G ( rQE,rQE, ω )on frequency determined by various multipolar plasmon modes excited by the QE, resulting in several peaks in the spectra of Im{G ( rQE,rQE, ω ) } (Figure 3.3a). The three lower-frequency sharp resonances are associated with the dipolar (DP, ℓ = 1), quadrupolar (QP, ℓ = 2), and octupolar (OP, ℓ = 3) plasmons of the spherical MNP. Their frequencies ωℓ are given by the poles of Eq. (3.9) , Reϵ(ωℓ) + ℓ+1 ℓ = 0. As 84 3.2. Results and discussion Figure 3.3:(a) Classical LRA result of the (a) imaginary part ( Im{G ( rQE,rQE, ω ) } ), and (b) real part ( Re{G ( rQE,rQE, ω ) } ) of the self-interaction Green’s function G ( rQE,rQE, ω )obtained for a point-like QE placed in front of an individual spherical MNP of radius a = 65 . 83 a0 ( ≈ 3 . 5 nm ). Results are shown as a function of the frequency ω of the oscillating QE and the surface-to-emitter distance, d . Panels (c) and (d) correspond to the results obtained with TDDFT simulations. In (a) and (b), the upper and lower range of values in the color bar denote saturation. we consider a metal described with a Drude dielectric function, this results in ωℓ=ωprℓ 2ℓ+ 1.(3.11) From Eq. (3.11) it follows that the frequencies of the DP, QP and OP are respectively ωDP ≈ 3 . 4 eV , ωQP ≈ 3 . 7 eV , and ωOP ≈ 3 . 85 eV . The highfrequency broad peak at ωPSM ∼ 4 eV (i.e., close to the surface plasmon frequency ωSP = ωp/√2≈ 4 . 16 eV ) corresponds to the so-called pseudomode [ 143 ], which is composed by a pilling up of several overlapping high-order plasmonic modes (ℓ= 4,5,6, . . . ) with closely-spaced resonant frequencies. 85 Chapter 3. Quantum surface effects in the electromagnetic coupling between quantum emitters and metallic nanoparticles At small d≈ 0 . 53 − 1 nm , Im{G ( rQE,rQE, ω ) } , as calculated within the classical LRA, is dominated by the pseudomode excitation. As the distance between the QE and the MNP increases, Im{G ( rQE,rQE, ω ) } decreases, and the relative contribution of different plasmon modes changes in favor of the low ℓ resonances (this behavior can be seen more clearly in Figure 3.4). Thus, for large d∼ 1 . 6 − 1 . 8 nm , the values of Im{G ( rQE,rQE, ω ) } attained within the pseudomode frequency range become comparable to those at the sharp DP and QP resonances. The faster decrease of the resonances associated with high-order plasmon modes with increasing d can be inferred from Eq. (3.9) , where R = a + d . We also note that the value of Im{G ( rQE,rQE, ω ) } obtained for a QE resonant with the pseudomode at ωPSM = 4 . 05 eV and located at d = 0 . 58 nm corresponds to a Purcell factor FP≈ 5 . 2 × 10 6 [Eq. (3.2) ]. This very large value is explained by the small volume of the MNP (and the resulting strong field localization). In contrast to the classical results, Im{G ( rQE,rQE, ω ) } as calculated with TDDFT (Figure 3.3c) mainly reveals a single broad feature [ 345 ] for the range of distances considered in this chapter. At small separation d≈ 0 . 53 − 1 nm , the maximum value of Im{G ( rQE,rQE, ω ) } is reached within the frequency interval ω∼ 3 . 6 − 3 . 7 eV , i.e., it is redshifted with respect to the classical pseudomode peak. As d increases, the resonant feature slightly shifts to lower frequencies. Moreover, the overall profile somewhat sharpens, albeit, in sheer contrast with the classical theory, the contributions of different plasmon modes remain spectrally broader and are barely resolved. Consistent with the strong broadening of the plasmon resonances due to quantum effects, the TDDFT results show smaller values of Im{G ( rQE,rQE, ω ) } at resonance, and thus lower QE decay rates, as compared to the classical LRA prediction. For example, for a QE placed at a distance d = 0 . 58 nm , the resonant Purcell factor FP calculated within TDDFT is FP≈ 1 . 5 × 10 6 . This is more than three times smaller than the maximum LRA value. On the other hand, the broadening of the spectra leads to a larger off-resonant Im{G ( rQE,rQE, ω ) } obtained with TDDFT as compared to classical LRA predictions. We next compare the classical LRA and quantum TDDFT results for the real part of the self-interaction Green’s function, Re{G ( rQE,rQE, ω ) } , which determines the Lamb shift ∆ωQE of the QE transition frequency [Eq. (3.1b)]. As depicted in Figure 3.3b, and consistent with the results obtained for the imaginary part of the Green’s function (Figure 3.3a), the frequency dependence of Re{G ( rQE,rQE, ω ) } obtained from classical LRA calculations features a rich resonance profile. For an individual plasmonic mode, the Kramers-Kronig relations would lead to a sign change of Re{G ( rQE,rQE, ω ) } at the resonance frequency. In the full calculations, Re{G ( rQE,rQE, ω ) } does not show the sign change at resonance for low ℓ modes and small distances, because of the off-resonant contribution associated with neighboring plasmon modes with larger ℓ . It is only at the pseudomode frequency that the contribution of the nearly degenerate resonances leads to a change of sign of Re{G ( rQE,rQE, ω ) } from positive values at frequencies below ωPSM ∼ 4 eV to negative values above this frequency. When d increases, the contribution from off-resonant neighboring modes is reduced so that, in addition to the pseudomode 86 3.2. Results and discussion resonance, the sign change of Re{G ( rQE,rQE, ω ) } can be observed at the DP and QP resonances. This appears particularly clear in Figure 3.4 discussed below, where we show the frequency dependence of the Green’s function calculated for a set of fixed separations, d, between the QE and the MNP surface. Similarly to the results obtained for the imaginary part, the TDDFT calculations in Figure 3.3d show smaller absolute values of Re{G ( rQE,rQE, ω ) } (Lamb shift) and a broader structure at the resonant plasmon frequency ω∼ 3 . 3 − 3 . 7 eV as compared to the classical LRA. This holds for the entire distance d range considered in this chapter. Notably, a single broad resonance is appreciated in the TDDFT results, and Re{G ( rQE,rQE, ω ) } changes its sign from positive to negative only at the resonant frequency ω∼ 3 . 3 − 3 . 7 eV , i.e. at lower frequency than within the classical LRA. As a consequence, for QEs with transition frequencies within the range of ωQE ∼ 3 . 7 − 4 eV , each model predicts a photonic Lamb shift ∆ ωQE ∝Re{G ( rQE,rQE, ω = ωQE ) } of opposite sign [Eq. (3.1b) ]. For example, according to LRA, a QE located at d = 0 . 58 nm and characterized by a transition dipole moment µQE = 0 . 1 enm and resonant frequency ωQE = 4 eV experiences a redshift of ∆ ωQE ≈ − 130 meV . In contrast, TDDFT predicts a blueshift of ∆ωQE ≈34 meV under the same conditions. 3.2.2 Quantum TDDFT vs. semiclassical nondispersive SRF After identifying the main quantum-mechanical effects in Figure 3.3, we can use the nondispersive SRF to dissect the role of nonlocality and the spill-out of the induced charges that can be behind the differences between the classical LRA and TDDFT results of G ( rQE,rQE, ω )discussed above. To this end, we compare in Figure 3.4 the real part (upper panels a-e) and the imaginary part (lower panels f-j) of G ( rQE,rQE, ω )as calculated using the three different approaches (TDDFT, LRA, and SRF). Results are shown as a function of the oscillation frequency of the QE, for selected values of the surface-to-emitter distance d . The solid and dashed blue lines show the reference TDDFT results and those obtained using the nondispersive SRF, respectively. The classical LRA results are plotted by gray-dotted lines. The overall good agreement between the TDDFT and the nondispersive SRF in Figure 3.4 establishes the validity of the latter and allows us to use the framework of the SRF to analyze the role of the quantum phenomena manifested in G ( rQE,rQE, ω ), as we discuss below. We first focus on the results at relatively large distance d = 1 . 38 − 2 . 22 nm (panels a-c and f-h), where the agreement between the nondispersive SRF and TDDFT is particularly good. The semiclassical SRF accurately reproduces the TDDFT results of the spectral position and resonance profile of G ( rQE,rQE, ω ), thus correctly accounting for the redshift and larger broadening of the peaks as compared to the classical LRA. On the other hand, for distances below d≤ 0 . 95 nm (panels d-e and i-j), the nondispersive SRF results are redshifted with respect to those of TDDFT, i.e., the semiclassical model based on the long-wavelength approximation of d⊥ ( ω )overestimates the redshift of the plasmonic modes from 87 Chapter 3. Quantum surface effects in the electromagnetic coupling between quantum emitters and metallic nanoparticles a cut of Figure 3.6. Although there are quantitative discrepancies between the TDDFT and the SRF results for narrow gaps D∼ 1 − 1 . 5 nm (panels d,e,i,j in Figure 3.7), the overall good agreement between the two approaches indicates that the spill-out of the induced charges and surface-enabled Landau damping, already discussed in the individual MNP results, are also the main quantum mechanisms influencing G ( rQE,rQE, ω )in the dimer configuration. Moreover, these results corroborate the validity of the long-wavelength limit implementation of the SRF to adequately describe the electromagnetic interaction between a QE and a plasmonic gap nanostructure for situations where the gap separation is larger than D∼ 1 . 5 nm (panels a,b,c,f,g,h in Figure 3.7). However, similar to the case of the individual MNP in Subsections 3.2.2 and 3.2.3, for smaller gap separations in a dimer, D < 1 . 5 nm , the limitations of the approximation used to obtain the Feibelman parameter d⊥(ω)(neglecting the nonlocality of optical response in the direction parallel to the metal surface) affects the accuracy of the results, since in such a case, high-order plasmonic ℓ -modes are also relevant in the response of the system. 3.3 Summary In this chapter, we study the influence of quantum phenomena on the electromagnetic interaction between a point-like quantum emitter (QE) and canonical metallic nanostructures. We focus on the study of the self-interaction dyadic Green’s function ˆ G ( rQE,rQE, ω )obtained for an individual spherical MNP and a dimer comprising two identical spherical MNPs, with the QE oriented perpendicular to the metal surfaces. In the case of the dimer, the QE is located in the middle of the gap. We consider sufficiently large QE–MNP separations so that charge-transfer processes related to electron tunneling do not play a role. We first calculate ˆ G ( rQE,rQE, ω )in the presence of an individual sodium MNP using time-dependent density functional theory (TDDFT), and then employ analytical expressions derived from a semiclassical model (referred to as the nondispersive SRF) in order to identify the origin of the quantum effects that influence the QE–MNP coupling. This nondispersive SRF incorporates surface quantum-response corrections by means of the Feibelman d⊥ ( ω )parameter obtained under the long-wavelength approximation. The overall good agreement between TDDFT and the nondispersive SRF for both the individual and the dimer configurations confirms that surface-enabled Landau damping and the spill-out of the induced electron density drastically affect the electromagnetic QE–MNPs interaction. These mechanisms explain why the resonances of ˆ G ( rQE,rQE, ω ) obtained from TDDFT are redshifted and broader as compared to those obtained from a classical calculation using the local-response approximation (LRA) of the optical response of the metals. We find that these quantum effects become more significant with increasing order ℓ of the plasmonic resonance. The analysis of the TDDFT calculations indicates that the higher the value of ℓ , the larger the broadening κℓ produced by 94 3.3. Summary surface-enabled Landau damping as well as the redshift produced by the spill-out of the dynamical screening charges. Thus, the quantum phenomena explored in this chapter show a considerable influence in the optical response for very small distance between a QE and a metallic surface, and when the QE is coupled to high-order plasmonic modes of the nanostructure. We also find that the nondispersive SRF calculations based on the longwavelength limit of the Feibelman parameter d⊥ ( ω )describe more accurately the Green’s function of a dipolar emitter when considering QEs coupled to loworder plasmonic modes. However, these calculations are not accurate when the contribution from high-order modes ( ℓ≳ 5) of the nanostructure is large, as occurs when the distance d between the QE and the MNP is very small, d∼ 0 . 6 nm . As a consequence, in the dimer configuration the nondispersive SRF calculations underestimate the broadening of the pseudomode and become inaccurate for gap separations of the order of D∼ 1 − 1 . 5 nm . The shortcomings of the nondispersive SRF are due to the limitations to address the parallel nonlocality of the response when using a d⊥ ( ω )parameter where the dependence on the wavenumber k∥ parallel to the metal surface (or, equivalently, the angular momentum ℓ for spherical MNPs) is neglected. For more accurate results, it is necessary to go beyond the longwavelength limit of d⊥ ( ω )to properly account for the nonlocality of the surface response in the direction parallel to the surface, as this impacts the results in situations where the QE–MNP distance is small (high-order ℓ modes involved). The extension of the SRF to a dispersive model that takes into account the k∥-dependence of d⊥is proposed in Chapter 4. This chapter thus provides a fundamental description of the quantum phenomena influencing the electromagnetic interaction between a QE and plasmonic nanostructures for surface-to-emitter distances as small as ≈ 0 . 5 nm . For even smaller separation, charge-transfer processes between the QE and the MNPs can influence the optoelectronic response of the system, so that a many-body treatment based on TDDFT of both the QE and the MNPs is necessarily required to naturally account for any quantum effect in those extreme situations, including also chargetransfer processes. We adopt this strategy in Chapter 5. 95 Chapter 4 DISPERSIVE SURFACE-RESPONSE FORMALISM TO ADDRESS OPTICAL NONLOCALITY IN SITUATIONS OF EXTREME PLASMONIC FIELD CONFINEMENT As shown in Chapter 3, the surface-response formalism (SRF) based on the Feibelman parameter d⊥ allows us to capture quantum effects such as the spill in/out of the induced charges, surface-enabled Landau damping, and nonlocal dynamical screening in a computationally simple manner. To the best of our knowledge, d⊥ has been so far computed using the long-wavelength approximation [ 98 ], which consists in neglecting the nonlocality of the optical response in the direction parallel to the metal–dielectric interface. The nondispersive Feibelman parameter is then a function of the excitation frequency, d⊥≡d⊥ ( ω ), and it does not depend on the wavenumber k∥ parallel to the metal surface. Considering the long-wavelength limit reduces the computational effort to obtain d⊥ ( ω )from quantum calculations. Moreover, it simplifies the implementation of the SRF in existing numerical tools that solve Maxwell´s equations in Nanophotonics as employed in a number of recent studies [97,124,129,296,299,300]. Using the long-wavelength limit of d⊥ (i.e., considering k∥∼ 0) is a reasonable approximation when the nonlocality of the optical response in the direction perpendicular to the surface dominates, i.e., when the characteristic scale ∆ s of the optical field variation along the surface is large as compared to that of the surface-response correction, ∆ s>> d⊥ . This is the case of e.g. typical individual metallic nanoparticles (MNPs) subjected to plane-wave illumination. However, in 97 Chapter 4. Dispersive surface-response formalism to address optical nonlocality in situations of extreme plasmonic field confinement Chapter 3we demonstrate that using the long-wavelength approximation of d⊥ ( ω ) within the nondispersive SRF is not accurate when describing situations involving multipole plasmon modes characterized by localized surface charges that rapidly vary along the MNPs surface. This is analogous to exciting plasmons with large transverse wavenumber k∥ , and therefore requires to go beyond the long-wavelength limit of d⊥. In this chapter, we demonstrate that accounting explicitly for the nonlocality of the optical response in the direction parallel to the surface lifts the shortcomings of the nondispersive SRF. Using a dispersive Feibelman parameter d⊥≡d⊥ ( ω, k∥ ) that is a function of ω and k∥ allows for correctly predicting the optical response of plasmonic structures in extreme situations where plasmon-induced charges characterized by high k∥ can be excited, such as in structures with small radius or small gaps, for instance. We first study in Subsection 4.2.1 the energy dispersion of localized multipolar plasmon resonances sustained by a cylindrical metallic nanowire using TDDFT, and show that the energies of multipolar plasmon modes in such a nanowire are governed by a universal parameter m/Rc equivalent to the wavenumber k∥ of surface plasmons at planar interfaces (here m is the magnetic quantum number, and Rc is the nanowire radius). Consistent with this statement, we demonstrate that the energy dispersion of localized multipolar plasmon modes relative to m/Rc in cylindrical nanowires follows the previously-studied k∥ -dispersion of nonretarded surface plasmons at planar metal surfaces [ 62 , 306 , 346 ]. Thus, TDDFT calculations of the cylindrical nanowire are used in Subsection 4.2.2 to obtain the Feibelman parameter d⊥≡d⊥ ( ω, k∥ = m/Rc ), which reveals a strong dependence on k∥ . We show in Subsection 4.2.3 that the dispersive ( k∥ -dependent) SRF accurately reproduces the energy dispersion of plasmon resonances obtained from TDDFT for the cylindrical nanowire. Furthermore, we demonstrate that the same set of d⊥ ( ω, k∥ )can be used to address optical nonlocality in the parallel direction in metallic nanostructures of different shapes. To this end, we apply the dispersive SRF to small spherical MNPs as well as to nanometer-gap spherical dimers coupled to quantum emitters (QEs), obtaining a good agreement between the dispersive SRF and TDDFT results. This chapter thus provides a significant advance toward the implementation of a SRF that adequately accounts for quantum effects in the optical response of plasmonic systems exhibiting extreme optical nonlocality. 4.1 System and methods Three different plasmonic nanostructures are considered in this chapter to show the generality of the dispersive SRF: (i) an infinite cylindrical nanowire of radius Rc extended along the z -axis, as described in Subsection 2.3.3, which is used to calculate d⊥ ( ω, k∥ ), (ii) an individual spherical MNP of radius a , as described in Subsection 2.3.1, and (iii) a dimer consisting of two identical spherical MNPs with a QE located in the middle of the gap. The dimer considered here is the same as the one considered in Chapter 3. All the nanostructures considered in this chapter are described within the jellium model of free-electron metals (see 98 4.1. System and methods Figure 4.1: Sketch of the systems studied in this chapter. (a) Cylindrical Na nanowire with radius Rc , infinite along the z -axis. ρ is the radial coordinate, and φ the azimuthal angle. The system possesses rotational symmetry with respect to the z -axis. We consider three different values of the radius: Rc = 75 a0, 100 a0,and 150 a0 .(b) Spherical Na nanoparticle with radius a = 65 . 83 a0 .(c) A point-like quantum emitter (QE) with a dipole moment pd placed at the middle of the gap formed by two identical spherical Na nanoparticles. The QE is represented by the green arrow. The gap separation distance is denoted as D . In all the structures, Na is characterized by a Wigner–Seitz radius rs= 4 a0. Subsection 2.1.2) using a Wigner–Seitz radius rs = 4 a0 that corresponds to sodium. The classical bulk plasma frequency is therefore ωp = q3 r3 s = 5 . 89 eV, the surface plasmon frequency ωSP = ωp/√2 = 4 . 16 eV, and the localized dipolar plasmon (DP) frequency ωDP = ωp/√3 = 3 . 4eV. We choose this material because it allows for a direct comparison between our results and those obtained in Chapter 3 using the long-wavelength approximation. The schematic representation of the nanostructures studied in this chapter is depicted in Figure 4.1. The optical response of individual spherical MNPs and dimers has been studied in detail in previous chapters, and thus here we focus on the numerical procedure used to obtain the energy dispersion of localized multipolar plasmon resonances sustained by the cylindrical nanowire. We address cylindrical nanowires with large radii Rc (as compared to that of spherical MNPs) within the range Rc = 75 − 150 a0 ( ≈ 4 − 8 nm )characterized by well-developed multipolar plasmons, thus allowing us to span a large range of periods 2 πRc/m of the spatial variation of plasmon-induced charges along the nanowire surface. Here, owing to the cylindrical symmetry of the system, the multipolar plasmon modes can be characterized by the magnetic number m related to the eimφ dependence on the azimuthal angle φ of the potentials, fields and induced charges. In the frequency domain, the potential Vind ( ρ, φ, ω )induced at ρ > Rc in response to an external excitation Vext(x, y, t)can be expressed as 99 Chapter 4. Dispersive surface-response formalism to address optical nonlocality in situations of extreme plasmonic field confinement Vind(ρ, φ, ω)=2 mmax X m=1 Rc ρm Qm(ω)cos(mφ),(4.1) where Qm ( ω )is the frequency-resolved multipole moment of order m per unit length along the z -axis. Our aim is to determine the frequency and width of these multipolar resonances. With this purpose, we use the Kohn–Sham scheme (KS) of TDDFT introduced in Section 2.2. We employ the same numerical implementation as described in Subsection 2.3.3 but, instead of considering a x -polarized plane-wave excitation given by Eq. (2.76) , here we consider the following impulsive external potential Vext(x, y, t): Vext(x, y, t) = ξ δ(t) mmax X m=1 ρ Rcm cos(mφ),(4.2) where x = ρcosφ and y = ρsinφ , δ ( t )is the Dirac delta function, the amplitude ξ is sufficiently weak to ensure a linear response, and mmax = 30 is the highest multipole order considered in the simulation. In Eq. (4.2) , we write explicitly the dependence of Vext ( x, y, t )on the spatial variables ( x, y )to stress that we use Cartesian coordinates as described in Subsection 2.3.3. Due to the excitation used here, multipolar plasmon excitations are localized at the ( x, y )-plane, and do not propagate along the z -axis. The term localized referring to plasmons excited at the nanowire is often omitted but implicitly assumed. Because of the cylindrical symmetry of the system, m is a good quantum number, i.e., an external potential with angular dependence cos ( mφ )excites localized multipolar plasmon modes at the nanowire characterized by induced charges with the same angular dependence cos(mφ). Using the real-time ALDA-TDDFT methodology in Cartesian coordinates as introduced in Subsection 2.3.3, we calculate the time-dependent multipole moment Qm ( t )induced at the nanostructure per unit length along the z -axis, defined as 25 Qm(t) = −1 mZZ dx dy ρ Rcm cos(mφ)δn(x, y, t),(4.3) where the electron density δn ( x, y, t )induced by Vext ( x, y, t )in Eq. (4.2) is given by [Eq. (2.77)] δn(x, y, t) = n(x, y, t)−n(x, y, t = 0).(4.4) The frequency-resolved spectrum of the multipole moment Qm ( ω )is finally obtained from the time-to-frequency Fourier transform, Qm(ω) = Zdt Qm(t)e(iω−η/2)t,(4.5) 25 The 1 /(Rc)m factor of the definition of the cylindrical multipole moment in Eq. (4.3) (and consistently the factor ( Rc ) m in Eq. (4.1) ) simplifies the comparison between the results of the calculations performed for nanowires of different radius Rc. 100 4.2. Results and discussion where an attenuation factor η = 0 . 15 eV is used [Eq. (2.56) ] to mimic dissipation processes beyond the reach of the ALDA-TDDFT scheme adopted here (see details in Section 2.2). It is worth mentioning that, for mmax = 1,Eq. (4.2) expresses the nonretarded potential corresponding to a plane-wave illumination polarized along the x -axis. Such illumination is typically used in linear-response TDDFT calculations to obtain the dipolar polarizability of the system, as described in Subsection 2.3.3. 4.2 Results and discussion 4.2.1 TDDFT study of the energy dispersion of multipolar plasmon resonances in a metallic nanowire Figure 4.2 shows the intensity spectrum of the multipole moments, |Qm ( ω ) |2 , obtained for a cylindrical nanowire of radius Rc = 150 a0 (panel a) and Rc = 100 a0 (panel b) in response to the external potential given by Eq. (4.2) . The results for different values of m are plotted, ranging from m = 1 (top) to m = mmax (bottom). For the largest nanowire mmax = 30 is considered, whereas for the smallest one mmax = 23. The general behavior of multipolar plasmon resonances |Qm ( ω ) |2 is independent of the size of the nanostructure: first, for a given m , a well-defined resonance centered at a frequency ωm is obtained associated with the excitation of the multipolar plasmon mode of order m . The width of the multipolar plasmon resonance for the two nanowires increases with increasing m because of the enhancement of surface-enabled Landau damping, where the plasmon decays into electron–hole pair excitations at the surface region, as discussed in Chapter 3 for spherical MNPs for increasing ℓ . Moreover, regardless of the specific value of the radius Rc , ωm first redshifts with increasing m and, after reaching a minimum, it continuously blueshifts. However, the multipolar plasmon resonances of a given order m are broader for the smaller nanowire, and do not emerge at the same frequency ωm for the two sizes. For example, the dipolar plasmon resonance ( m = 1) emerges at ω1∼ 4 . 1 eV for Rc = 150 a0 (extremely close to the classical long-wavelength limit of the surface plasmon frequency, ωSP = 4 . 16 eV ), while for Rc = 100 a0 it appears at ω1∼ 4 . 06 eV . Moreover, for the largest nanowire the minimum value of ωm is reached at m = 10 ( ω10 ∼ 3 . 8 eV ), whereas for the smallest nanowire the minimum is obtained at m = 8. The differences in the results obtained for the two nanowires are particularly apparent when comparing ωm for m= 20: in this case, there is a mismatch in ωmof the order of ∼0.5eV. The results shown in Figure 4.2a,b are closely related to the dispersion relation of surface plasmon resonances supported by planar metal–vacuum interfaces as a function of the parallel wavenumber k∥ , as we discuss below. Indeed, the induced fields and surface charge densities of localized multipolar plasmons considered in this chapter have a dependence of exp ( imφ )on the azimuthal angle φ . By introducing the coordinate r∥ along the surface of the nanowire cross-section, r∥ = Rcφ , the angular dependence transforms into exp ( imφ ) →exp ( im Rcr∥ ). 101 Chapter 4. Dispersive surface-response formalism to address optical nonlocality in situations of extreme plasmonic field confinement Figure 4.2:(a) Intensity spectrum of the multipole moment |Qm ( ω ) |2 induced at an infinitely long cylindrical Na nanowire ( rs = 4 a0 ) of radius Rc = 150 a0 ( ≈ 8 nm ). Results are shown as a function of the frequency of the external excitation, ω , for different values of the magnetic quantum number m , ranging from m = 1 (top) to m = 30 (bottom), as indicated in the insets. All |Qm ( ω ) |2 are normalized to their corresponding maximum value. (b) Same as in (a) but for a smaller nanowire with radius Rc = 100 a0 ( ≈ 5 . 3 nm )and magnetic number within the range m= 1 −23. Results are displaced in the vertical axis for clarity. For Rc→ ∞ , the cylindrical geometry tends to the planar-surface geometry with r∥ being the coordinate parallel to the surface, so that setting m Rc→k∥ recovers the standard dependence exp ( ik∥r∥ )of a surface plasmon propagating along the surface [see Eq. (1.38) ]. Thus, we interpret m Rc in the following as an “effective” wavenumber k∥ by considering that localized multipolar plasmons in the nanowire correspond to confined surface plasmons with a quantized wavelength λ∥ = 2 πRc/m [ 347 ]. To support this correspondence, we compare in Figure 4.3 the intensity spectra |Qm ( ω ) |2 obtained for nanowires with radius Rc = 75 a0 ( dahsed lines ) , 100 a0 ( dotted lines ) ,and 150 a0 ( solid lines ). Results are shown for selected values of m such that k∥ = m/Rc = 0 . 013 , 0 . 02 , 0 . 04 , 0 . 067 ,and 0 . 1 a−1 0 . Whenever the ratio m/Rc is fixed, the intensity spectrum of the multipole moment |Qm ( ω ) |2 obtained for nanowires of different size shows near perfect match with each other, thus confirming that an effective wavenumber k∥=m/Rcdetermines the optical response of the system. 102 4.2. Results and discussion Figure 4.3: Intensity spectrum of the multipole moments |Qm ( ω ) |2 obtained for cylindrical nanowires of different radii Rc and selected values of the magnetic number m . Solid lines: Rc = 150 a0 . Dotted lines: Rc = 100 a0 . Dashed lines: Rc = 75 a0 . The selected values of m are such that k∥ = m/Rc = 0 . 013 a−1 0 (red), k∥ = 0 . 02 a−1 0 (blue), k∥ = 0 . 04 a−1 0 (brown), k∥ = 0 . 067 a−1 0 (orange), and k∥ = 0 . 1 a−1 0 (purple) for the three values of Rc . The spectra |Qm ( ω ) |2 corresponding to a specific wavenumber k∥ = m/Rc are normalized to the maximum value obtained for the case of Rc= 150 a0. Results are displaced in the vertical axis for clarity. 4.2.2 Calculation of the dispersive Feibelman parameter d⊥(ω, k∥) The results presented in the previous section allow us to calculate d⊥≡d⊥ ( ω, k∥ )as a function of both the excitation frequency ω and the effective wavenumber parallel to the surface, k∥ = m/Rc , using the cylindrical geometry. In this subsection, we first obtain the expression of the Feibelman parameter d⊥ for a cylindrical metallic nanowire within the SRF, and then discuss the TDDFT results of d⊥≡d⊥ ( ω, k∥ ) calculated for the present system following the methodology employed in the previous subsection. Expression of the Feibelman parameter d⊥for a cylindrical nanowire To obtain the expression of the Feibelman parameter d⊥ for a cylindrical nanowire infinite along the z -axis, we focus on the SRF solution of the electrostatic potential ϕind ( ρ, φ, ω )induced at the nanowire (in the nonretarded approximation). Due to the translational invariance of the system with respect to the z -axis, ϕind only depends upon the spatial variables ( ρ, φ ). The induced potential ϕind ( ρ, φ, ω )can thus be expressed as: ϕind(ρ, φ, ω) = m=∞ X m=−∞ ϕm(ρ, ω)eimφ,(4.6) where formally the sum extends from m=−∞ to m=∞. The radial part of the induced potential, ϕm ( ρ, ω ), is given by the solution of 103 Chapter 4. Dispersive surface-response formalism to address optical nonlocality in situations of extreme plasmonic field confinement D=2.33 nm D=2.33 nm D=1.06 nm D=1.06 nm Figure 4.6: Comparison between the results obtained using TDDFT (solid lines), dispersive SRF (long-dashed lines), and nondispersive SRF (short-dashed lines). (a) Imaginary part of the first ten multipolar polarizabilities αℓ ( ω )( ℓ = 1 − 10). The left-hand side panel presents the comparison between TDDFT and nondispersive SRF results, whereas the right-hand side panel presents the comparison between TDDFT and dispersive SRF results. Each spectra is normalized to the corresponding maximum value obtained within TDDFT for each value of ℓ . The spectra corresponding to different ℓ are vertically displaced for visibility. The TDDFT results are represented by solid lines with hatched area. (b,c) Lamb shift ∆ ωQE (left-hand side panels) and Purcell factor FP (right-hand side panels) obtained within the three methods for a point-dipole quantum emitter (QE) at the center of a spherical MNP dimer of radius a = 65 . 83 a0 ( ≈ 3 . 5 nm ). The dipole is oriented along the dimer axis, and its transition dipole moment is µ = 0 . 1 enm (with ethe electron charge). In (b), the gap separation is D= 2.33 nm. In (c),D= 1.06 nm. 110 4.2. Results and discussion plasmon mode appears in each spectrum of Im{αℓ} , which allows for the discussion on their energies and widths. The TDDFT results of Im{αℓ ( ω ) } in Figure 4.6a show that the multipole plasmon resonances of order ℓ continuously blueshift with increasing ℓ in the considered range ℓ = 1 − 10. The resonance broadens as ℓ increases due to the enhancement of surface-enabled Landau damping [ 64 , 66 , 194 , 259 – 261 ]. As already discussed in Chapter 3, the nondispersive SRF accurately reproduces the TDDFT data of Im{αℓ ( ω ) } for low values of ℓ∼ 1 − 4, but fails to make correct predictions for ℓ≥ 5. Indeed, for these large values of effective k∥ = ℓ/a , the plasmonic resonances within the nondispersive SRF start to redshift with increasing ℓ in contrast to the continuous blueshift obtained from TDDFT calculations. Thus, the nondispersive SRF predicts multipole plasmon frequencies that deviate significantly from the TDDFT values. In sheer contrast, by accounting for the dependence of the Feibelman parameter d⊥ on k∥ , the dispersive SRF correctly captures the energy blueshift and broadening of multipole plasmon resonances in Im{αℓ ( ω ) } (see right-hand side panel in Figure 4.6a). Although some quantitative differences emerge for large multipole order ℓ = 7 − 10, one can observe an overall good agreement between TDDFT and the dispersive SRF results over the entire range of ℓ values considered here. Thus, the dispersive SRF is useful to describe localized multipole plasmon resonances of large order ℓsustained by small MNPs. Finally, we address another canonical plasmonic system: a dimer of spherical MNPs. Specifically, we study the case of a point-dipole quantum emitter (QE) located at the center of the gap formed by two identical spherical MNPs with radius a = 65 . 83 a0 ( ≈ 3 . 5 nm ), as sketched in Figure 4.1c. This system is identical to the one considered in Chapter 3. The QE is oriented along the axis of the MNP dimer (the z -axis). The gap separation distance, D , is in the nanometer scale, and thus nonlocality strongly influences the optical response of the system, as shown in Chapter 3. We focus on the enhancement of the QE total decay rate given by the Purcell factor FP and the change of resonant frequency ∆ ωQE (Lamb shift) due to the self-interaction of the QE with the MNP dimer. The Lamb shift is calculated considering a transition dipole moment µ = 0 . 1 enm (with e the electron charge). The TDDFT and SRF results are obtained within the nonretarded approximation following the procedure described in Chapter 3, where for the dispersive SRF we use the Feibelman parameter d⊥(ω, k∥)obtained in this chapter. Figure 4.6b shows the Lamb shift ∆ ωQE (left-hand side panel) and Purcell factor FP (right-hand side panel) obtained for a gap separation D = 2 . 33 nm , as calculated with the three models employed in this chapter (TDDFT, dispersive SRF, and nondispersive SRF). The three approximations show qualitatively good agreement, although the nondispersive SRF results slightly deviate from the TDDFT and dispersive SRF predictions. For this relatively large gap, the excitation of lowℓ multipole plasmon resonances dominates the response of the MNP dimer to the field created by the point-dipole QE (see Chapter 3), which validates the longwavelength approximation behind the nondispersive SRF results. Nonetheless, the results obtained within the dispersive SRF are more accurate when compared to TDDFT. 111 Chapter 4. Dispersive surface-response formalism to address optical nonlocality in situations of extreme plasmonic field confinement The better performance of the dispersive SRF to describe the electromagnetic QE–MNPs interaction is more evident when considering a smaller gap, which naturally involves larger values of k∥ in the response. Figure 4.6c shows the Lamb shift ∆ ωQE and Purcell factor FP for a gap separation D = 1 . 06 nm . In this situation, because of the higher spatial confinement of the induced charges at the metal surfaces across the gap, plasmon modes with large multipolar order ℓ become important. These largeℓ modes have overlapping resonant frequencies and thus contribute to a single broad peak (referred to as the pseudomode, see Chapter 3) at ω∼ 3 . 4 eV , as revealed by the TDDFT calculations. Since the nondispersive model does not accurately describe the energy of largeℓ multipolar modes for the individual MNP (Figure 4.6a), it also fails to predict the energy and the width of the plasmon pseudomode obtained within TDDFT for the dimer of small gap. Moreover, the nondispersive SRF strongly overestimates the Purcell factor and the Lamb shift close to the bonding dipolar plasmon (BDP) resonance at ω∼ 2 . 75 eV because of the contribution of highℓ multipolar modes near the BDP frequency within the nondispersive model (see Figure 4.6a). In contrast, the dispersive SRF provides accurate results even for this small gap separation, thus indicating that the dispersive SRF is well suited to correctly account for nonlocality in situations where plasmon-induced charges are characterized by a rapid variation in the direction parallel to the metal surface. 4.3 Summary In summary, in this chapter we have proposed a dispersive SRF that explicitly accounts for the dependence of the Feibelman parameter d⊥ on the wavenumber parallel to the metal surface, k∥ . Using TDDFT calculations as a reference, we have demonstrated that the dispersive SRF is much more accurate than the nondispersive SRF, usually implemented in the literature, in describing plasmonic systems characterized by extremely confined induced fields. The dispersive SRF proposed here thus overcomes the limitations of the nondispersive SRF identified in Chapter 3. Using the analogy between localized multipolar plasmons in infinite cylindrical nanowires of radius Rc and propagating surface plasmons at planar metal– vacuum interfaces, we have demonstrated that m/Rc can be interpreted as a wavenumber parallel to the surface, k∥ = m/Rc (here m is the magnetic quantum number). This study has allowed us to obtain the dispersive Feibelman parameter d⊥ ( ω, k∥ = m/Rc )using cylindrical nanowires, which is incorporated into the SRF to complete the description based on the long-wavelength value d⊥ ( ω, k∥ = 0), used in Chapter 3and in other recent works [124,129,296,299,300]. Supported by the examples of cylindrical and spherical metallic nanostructures, we have demonstrated that, in contrast to the nondispersive model, the dispersive SRF accurately describes the nonlocal optical response in extreme situations where the induced charges are characterized by a rapid variation in the direction parallel 112 4.3. Summary to the metal surface (large k∥ ). The results shown in this chapter thus contribute to the development of a theoretical model that captures quantum nonlocal effects in extreme situations, while keeping the numerical efficiency and easy implementation into the framework of classical electromagnetic theories [ 348 ]. We thus believe that the dispersive SRF proposed in this chapter can be useful to correctly account for optical nonlocality in nanostructured systems with extreme plasmonic field confinement, as it can be the case of metallic nanostructures interacting with fast electrons, MNPs coupled to QEs in close proximity, or MNPs ensembles with extremely narrow junctions. The dispersive SRF substantially improves the perfomance of the nondispersive SRF to describe nonlocality in the optical response of narrow junctions, however it still lacks the description of charge-transfer processes. To account for such effects by using semiclassical models, it would be necessary to further develop the present framework combining the SRF with e.g. a quantum-corrected model [83–87]. 113 Chapter 5 ELECTRONIC EXCITON–PLASMON COUPLING IN A NANOCAVITY BEYOND THE ELECTROMAGNETIC INTERACTION PICTURE The main mechanism that controls the interaction between quantum emitters (QEs) and metallic nanoparticles (MNPs) is the excitation of the QE exciton by the local electric field associated to the MNP plasmon. In previous chapters, the plasmonic response of MNPs is described either classically (Chapter 1) or by using TDDFT simulations that capture nonlocal and quantum surface effects (Chapter 3), while the exciton dynamics of the QE is modeled within the point-dipole approximation. The success of the methodologies used in previous chapters to explain the main features of the optical response in plasmonic nanocavities is due to the dominance of the electromagnetic interaction in the QE–MNP coupling for separations as small as one nanometer [ 145 , 349 ]. However, at even smaller separations between emitters and metal surfaces, of the order of Ångstroms, another quantum effect becomes important: electronic states localized at the QE and at the MNPs hybridize into "supermolecular" states which modify optical transitions, allowing for electron transfer between the QE and the MNP. Despite its importance [ 166 , 349 – 352 ], the effect of hybridization between the QE and the MNP electronic states as well as the corresponding electron-transfer processes remain largely unexplored in Nanophotonics, as the quantum theoretical treatment of the problem is challenging. It is only recently that such studies have become within the reach of theoretical efforts [ 107 – 109 , 323 , 329 ] enabling e.g. a better understanding of light emission in tunneling junctions [ 309 – 311 ]. Notably, it has been shown that a QE bridging two MNPs can trigger electron conductance 115 Chapter 5. Electronic exciton–plasmon coupling in a nanocavity beyond the electromagnetic interaction picture across subnanometric junctions, which strongly influences the optic and electronic (optoelectronic) response of the coupled system [107–109,329,353]. In this chapter, we apply a fully quantum many-body approach based on TDDFT to study the optoelectronic response and exciton dynamics in a QE–MNPs system where the QE is located at subnanometric separation from the metallic interfaces. In contrast with the methodologies employed in previous chapters, here we use a TDDFT treatment to describe the electronic structure of both the QE and the MNPs. We place particular emphasis on the role of electronic coupling and electron transfer between the QE and the MNPs to unveil the manifestation of these quantum effects in the optical response of the entire coupled system. Importantly, we demonstrate that the modification of the electronic structure of the hybrid QE–MNPs system as well as the broadening of the electronic states of the QE due to charge transfer lead to a breakdown of the classical electromagnetic description of plasmon–exciton interaction. We reveal important quantitative and qualitative differences between quantum TDDFT and classical LRA results of the linewidths and frequencies of the relevant optical modes. Moreover, we also observe the formation of a novel charge-transfer plasmon mode at low frequencies mediated by the emitter electronic structure. 5.1 System and methods We consider a QE interacting with a plasmonic dimer formed by two spherical MNPs. As sketched in Figure 5.1a, the QE placed in the middle of a plasmonic nanogap is illuminated by a plane wave polarized along the dimer axis ( z -axis). In this chapter, the gap separation D is varied to explore different regimes of electronic QE–MNPs coupling, ranging from electronically decoupled QE–MNPs (large D ) to electronically coupled ones (small D ). The calculation of the optical response is performed within the Kohn–Sham (KS) scheme of time-dependent density functional theory [ 70 – 72 , 189 ] (TDDFT) as introduced in Section 2.2, which successfully incorporates quantum phenomena such as many-body and single electron–hole pair excitations, electronic spill-out, nonlocal screening or electron tunneling in (sub)-nanometric metallic cavities [79,125,202,211,212,214,337,338,342]. The electronic structure of the MNPs is described within the jellium model of free-electron metals [ 194 , 200 ] introduced in Subsection 2.1.2, using a Wigner– Seitz radius of rs = 4 a0 that corresponds to sodium. Each MNP contains 638 conduction electrons (radius a = 34 . 4 a0≈ 1 . 8 nm ), and the Fermi level of the MNPs stands at EF = − 2 . 86 eV below the vacuum level. In contrast to the point-dipole approximation employed in previous chapters to model the QE (e.g. in Chapter 3), here we consider a “more realistic” QE that has a finite spatial extension. The electronic structure of the QE is described as a two-level system using a model potential VQE ( r )(see below). The optical response of the coupled QE–MNPs is addressed using the wave-packet propagation (WPP) method in cylindrical coordinate system, as introduced in Subsection 2.3.2. To excite the system, we apply an external potential Vext ( r, t ) = E0 ∆ t z δ ( t )[Eq. (2.67) ], which 116 5.1. System and methods 0 0.04 0.08 0.12 0.16 0=1 eV V0=3 eV V0=5 eV V Figure 5.1:(a) Sketch of the system studied in this chapter. A QE with a single optically-allowed HOMO–LUMO transition is located in the middle of a gap of size D formed by two spherical MNPs. The QE and each MNP contain 2 and 638 conduction electrons, respectively. (b) Effective one-electron potential Veff ( r )(top) and equilibrium electron density n0 ( r )(bottom) along the symmetry z -axis for the coupled QE–MNPs system with gap size D = 26 a0 . The HOMO and LUMO energy levels of the isolated QE are represented by red and green lines, respectively. The Fermi level EF = − 2 . 86 eV of the MNPs is shown by the black dashed line. (c) Absorption cross-section spectra σabs ( ω )of the isolated MNP dimer for D = 26 a0 (blue line) and D = 38 a0 (green). Dashed and solid lines correspond to the results obtained with classical LRA and TDDFT simulations, respectively. (d) Absorption cross-section spectra σabs ( ω )of the isolated QE for different values of the parameter V0 [Eq. (5.3) ] used to control optical and electronic properties of the QE. corresponds to plane-wave illumination polarized along the z-axis. Prior to studying the coupled QE–MNPs system, we summarize in Figure 5.1c the TDDFT results of the absorption spectra of the isolated MNP dimer of radius a = 34 . 4 a0≈ 1 . 8 nm considered in this chapter (solid lines). Gap separation distance of D = 38 a0 (green) and D = 26 a0 (blue) are considered here. The optical response of the MNP dimer is characterized by a bonding dipolar plasmon (BDP) resonance at ωBDP ∼ 3 eV . As expected from the results shown in Chapter 1 and Chapter 2, this BDP mode redshifts when reducing the gap separation because of the increased capacitive coupling between the two MNPs [ 354 ]. As a reference, in Figure 5.1c we also show the results from classical (nonretarded) LRA calculations (dashed lines) introduced in Section 1.1. The MNPs are described in this case with a Drude dielectric function [Eq. (1.10) ] using an “effective” plasma frequency ωp = 5 . 43 eV and intrinsic damping parameter γp = 0 . 15 eV . These parameters are chosen to provide reasonably good agreement between the TDDFT and classical 117 Chapter 5. Electronic exciton–plasmon coupling in a nanocavity beyond the electromagnetic interaction picture LRA results of the absorption spectrum of the individual MNP (not shown) [ 125 ]. Using this value of ωp = 5 . 43 eV (i.e., slightly smaller than the nominal bulk plasma frequency ωp = 5 . 89 eV of sodium for rs = 4 a0 ,Eq. (1.12) ) allows us to account for the redshift of the dipolar plasmon (DP) frequency of the small individual MNP because of the electron spill-out and dynamical screening introduced in Subection 2.3.1. In addition, the damping parameter γp = 0 . 15 eV used in LRA in this chapter accounts for all the decay channels of the plasmon excitation including the contribution of surface-enabled Landau damping. Further, in all classical LRA calculations we also introduce a gap scaling of ∆ = 3 . 4 a0 to (partially) account for the spill-out of the induced electron density with respect to the geometrical surface of the MNPs in the dimer configuration (correctly captured by the TDDFT simulations, as shown in Subsection 2.3.2). Introducing the gap scaling ∆=3 . 4 a0 is similar to considering the Feibelman parameter d⊥ ( ω )employed in Chapter 3and Chapter 4, as reported in ref. 58 . In this thesis we found that this procedure is valid to reproduce the redshift of the BDP of the MNP dimer, however it overstimates the strength of the higher-order plasmon modes such as the BQP. Thus, nonlocal and quantum effects are partially introduced in the classical LRA calculations of this chapter in an effective manner, which allows us to correctly reproduce the TDDFT spectra of the isolated MNP dimer. We next introduce in Subsection 5.1.1 the TDDFT description adopted in this chapter to model the QE, which allows for studying the effect of the electronic interaction between the QE and MNPs states. 5.1.1 Characterization of the model quantum emitter (QE) The electronic structure of the QE is described within the free-electron jellium model (Subsection 2.1.2), in a similar way as we model the MNPs. We consider a spherical QE of radius RQE = 5 a0 ( ≈ 0 . 26 nm )containing two valence electrons. The spin-restricted case is considered [ 355 ]. In the ground-state configuration of the QE, the total spin is zero, and the 2 electrons with opposite spins occupy the same Kohn–Sham (KS) valence orbital. From the charge neutrality condition, the positive background density n+ [Eq. (2.13) ] representing the atomic cores of the QE and spreading over its spatial extent satisfies 4 3πR3 QEn+= 2.(5.1) The occupied ( j = 1) and unoccupied ( j = 2 , 3 , 4 , . . . ) one-electron KS orbitals Ψ 0 QE,j ( r )of the QE and their energies ϵQE,j are obtained from the time-independent KS equation of DFT [Eq. (2.5)], ˆ H[n0 QE(r)]Ψ0 QE,j(r) = ϵQE,jΨ0 QE,j(r),(5.2) where ϵQE,j are the one-electron energy levels of the QE, and we use spherical coordinates as described in Subsection 2.3.1. The equilibrium electron density of the isolated QE, n0 QE ( r ), is given in this case by n0 QE ( r )=2 | Ψ 0 QE,1 ( r ) |2 , with the 118 5.1. System and methods factor 2 accounting for spin degeneracy. Controlling the values of the one-electron energy levels ϵQE,j of the QE states allows us to study different situations of the QE–MNPs coupling, where the QE exciton can be either in resonance or out of resonance with the main BDP resonance of the MNP dimer at ωBDP ∼ 3 eV (see Figure 5.1c). To this end, we introduce in the Hamiltonian ˆ H [ n0 QE ( r )] [Eq. (5.2) ] an additional attractive potential VQE ( r ) [ 216 ]. This attractive potential VQE ( r )can be thought of as a pseudopotential due to the atomic cores, it is localized in the spatial region of the QE, and it is given by VQE(r) = −V0e−4|r|2/R2 QE ,(5.3) The one-electron energy levels ϵQE,j of the QE states can be thus modified by changing the parameter V0. In this chapter, we refer to Ψ 0 QE,1 ( r )as the highest occupied molecular orbital (HOMO), which is a 1 s ( ℓ = 0 , m = 0) orbital with zero orbital momentum ( ℓ ) and magnetic quantum number ( m ). Thus, the electronic configuration of the QE is 1 s2 . The energy level of the HOMO is EHOMO = ϵQE,1 . For the values of V0 [Eq. (5.3) ] considered in this chapter, we find only three energy-degenerate unoccupied KS orbitals accessible for optical transitions from the ground state. These orbitals correspond to the 2 p -shell and are characterized by the orbital momentum ℓ = 1 and magnetic quantum numbers m = 0 ,± 1. With ψQE,2p ( r )the radial part of the KS orbital of the isolated QE, we can define Ψ0 QE,2(r) = Y0 1(θ, φ)ψQE,2p(r), Ψ0 QE,3(r) = Y−1 1(θ, φ)ψQE,2p(r), Ψ0 QE,4(r) = Y1 1(θ, φ)ψQE,2p(r), (5.4) where Ym ℓ ( θ, φ )are the spherical harmonics. Because of the symmetry of the system considered in this chapter (see Figure 5.1a), with an incident electromagnetic wave polarized along the z -axis, the electronic transitions preserve the magnetic quantum number m , and thus are effective between the ground-state 1 s KS orbital and the 2 p ( ℓ = 1 , m = 0) KS orbital Ψ 0 QE,2 ( r ). Thus, the optical absorption of the QE is determined by the 1 s→ 2 p transition. For the sake of simplicity, in this chapter we refer to Ψ0 QE,2(r)as the lowest unoccupied molecular orbital (LUMO). The free parameter V0 of the potential VQE ( r )[Eq. (5.3) ] is used to control the energy levels of the HOMO ( EHOMO = ϵQE,1 ) and the LUMO ( ELUMO = ϵQE,2 ) of the QE, as schematically depicted in Figure 5.1a by the red and green lines, respectively. As a consequence, the parameter V0 also determines the oscillator strength α0 QE and transition frequency ωQE of the QE exciton [Eq. (1.54) ] relevant in the optical response of the coupled QE–MNPs system. We show in Table 5.1 the energy levels EHOMO and ELUMO , as well as the oscillator strength α0 QE and transition frequency ωQE of the QE exciton, obtained for the three different values of the background potential V0 = 1 eV, 3 eV,and 5 eV [Eq. (5.3) ] considered in this chapter. The values of EHOMO and ELUMO are directly obtained from time-independent DFT calculations following the procedure described in 119 Chapter 5. Electronic exciton–plasmon coupling in a nanocavity beyond the electromagnetic interaction picture The classical LRA results of σabs ( ω )in Figure 5.5a show a splitting between the LR (blue dots) and the UR (green dots) with respect to the resonant frequency ωQE of the isolated QE already for large gap separation distance D = 40 a0 . This LR–UR splitting is a signature of the strong coupling between the QE exciton and the BDP resonance of the MNP dimer, as we further confirm by analyzing in Figure 5.5b the coupling strength g . The coupling strength g is obtained by fitting the classical LRA results of σabs ( ω )to the spectra obtained from the coupled harmonicoscillator model introduced in Section 1.4.3 [Eq. (1.60) ]. Figure 5.5b shows that, for the largest distance considered ( D = 40 a0 ), the criterion g > ( γMNP + γQE ) / 4 often used to identify strong coupling is satisfied. With decreasing D , the LR–UR frequency difference in Figure 5.5a strongly increases within the classical LRA model owing to the stronger electromagnetic coupling between the exciton and the plasmon. This is consistent with the increase of g , observed in Figure 5.5b, which satisfies the more restrictive criterion g > ( γMNP + γQE ) / 2for smaller gaps D≤ 32 a0 . Moreover, within the classical LRA framework, the higher-order plasmonic modes of the MNPs contribute to the electromagnetic interactions between the dimer and the QE for small gap separations. The effect of these higher-order modes, well documented for isolated dimer antennas [ 23 ], results here in an additional redshift of both the LR and UR branches, which explains why the UR branch appears at lower frequencies than the exciton frequency ωQE = 2 . 95 eV of the isolated QE. Results in Figure 5.5a,b thus demonstrate that the resonant QE–MNPs system would be in the strong-coupling regime according to the classical LRA description. The TDDFT results of σabs ( ω )in Figure 5.5c also show a splitting between the LR and the UR with respect to ωQE for large gap separation distance D = 40 a0 , although the strength of the LR is weaker than the one predicted by the classical model. This difference is a consequence of the nonlocality and finite-size effects (introduced in Section 2.3) that affect the optical response of the isolated MNP dimer, as confirmed by the analysis shown in Subsection 5.2.2 below. The general similarity between classical and TDDFT results for large separation occurs because there is no QE exciton quenching produced by electronic hybridization for such large separation. At those distances the electronic QE–MNPs coupling does not play a role. However, upon reducing the gap size D, the electronic orbitals of the QE hybridize with those of the MNPs. As shown in Figure 5.5d, the LUMO evolves into a broad structure reflecting the fast transfer of the excited electron between the LUMO of the QE and the conduction-band states of the MNPs quantized by the finite-size effect. This electronic interaction has an immediate consequence on the optical response of the strongly coupled QE–MNPs system since it hinders the energy transfer between the QE and the MNPs, thus attenuating the UR–LR splitting in exciton-plasmon polariton systems as well as producing a progressive merging of the LR and UR branches into a broad spectral feature when decreasing gap separation D . In this situation, strong electromagnetic coupling is, therefore, frustrated due to electronic QE–MNPs coupling. 126 5.2. Results and discussion 34 a0 38 a0 30 a0 26 a0 22 a0 18 a0 34 a0 38 a0 30 a0 26 a0 22 a0 18 a0 Figure 5.6: Absorption cross-section spectra σabs ( ω )of the studied QE–MNPs system, as obtained from the semiclassical model employed in this section. This semiclassical approach is based on the self-interaction Green’s function [Eq. (1.57) ] and considers that the QE is a classical point dipole, however ˆ G ( rQE,rQE, ω ), Eind ( r = rQE, ω ), ˆαMNP ( ω )and ˆαQE MNP ( ω )are obtained from TDDFT simulations of the isolated MNP dimer. Results are shown as a function of the frequency ω for gap size ranging from D = 16 a0 to D = 40 a0 in steps of 2 a0 . Panel (a) corresponds to the results obtained for ωQE = 2.58 eV, and panel (b) to ωQE = 2.95 eV. 5.2.2 Quantum finite-size effects vs. electronic QE–MNPs coupling In order to gain a better understanding of the role played by the electronic QE– MNPs coupling in the optical response, and to discard nonlocal and finite-size effects as studied in Chapter 3in connection with QE exciton quenching, we apply here a semiclassical 27 approach to the current QE–MNPs system. This semiclassical approach adopts the self-interaction Green’s function formalism (see Subsection 1.4.2), where the QE is introduced as a classical point dipole. The dipole moments induced at the MNPs and at the QE are then obtained from Eq. (1.57). However, the quantities ˆ G ( rQE,rQE, ω ), Eind ( r = rQE, ω ), ˆαMNP ( ω )and ˆαQE MNP ( ω ) are obtained from the TDDFT simulations of the MNP dimer, as described in Subsection 2.3.2 and Chapter 3. This semiclassical approach naturally includes finite-size effects on the response of the MNP dimer such as electron spill-out, nonlocality, surface-enabled Landau damping, and single electron–hole transitions. However, since the QE is introduced as a classical point dipole, the electronic coupling between the MNPs and the QE, as well as the actual electronic structure of the QE and the finite-size extension of its transition density are not accounted for. Thus, we expect that the differences between TDDFT and the semiclassical approach reveal the effect of electronic hybridization, only accounted for within the fully quantum TDDFT model. We show in Figure 5.6 the absorption cross-section σabs ( ω )obtained within the semiclassical model employed in this subsection for ωQE = 2 . 58 eV (panel a) 27 The semiclassical approach employed in this chapter should not be confused with the semiclassical SRF employed in Chapter 3and Chapter 4. 127 Chapter 5. Electronic exciton–plasmon coupling in a nanocavity beyond the electromagnetic interaction picture and ωQE = 2 . 95 eV (panel b). For large gap separation distances, D∼ 30 − 40 a0 , the semiclassical results shown in Figure 5.6 coincide with the results of TDDFT calculations displayed in Figure 5.2b ( ωQE = 2 . 58 eV ) and Figure 5.5c ( ωQE = 2 . 95 eV ). In particular, consistent with the discussion of the previous subsection, the LR calculated for the resonant case ωQE = 2 . 95 eV within both the semiclassical model and TDDFT is considerably weaker than the one predicted by the classical LRA approach (Figure 5.5a). Thus, the origin of the difference between the TDDFT and the classical LRA absorption spectra resides on quantum surface effects that are important for such small MNPs. Indeed, at large D there is no hybridization between the electronic states localized at the QE and at the MNPs. As already discussed, the electronic QE–MNPs coupling strongly affects the absorption spectra of the system for D below D∼ 26 a0 , which is now further corroborated from the comparison between the semiclassical results in Figure 5.6, the TDDFT results in Figure 5.2b and Figure 5.5c, and the classical results in Figure 5.2a and Figure 5.5a. First, TDDFT shows substantial broadening and reduction of the amplitude for the UR evolving from the BDP of the MNP dimer as compared to both semiclassical and classical results. These effects, not captured by the semiclassical model, are attributed to the charge-transfer processes between the MNPs. In our system, electron transport can occur at larger gap separations as compared to typical vacuum junctions (see Subsection 2.3.2) because it is assisted by photoexcited electron transfer through the LUMO of the QE [356]. On the other hand, the semiclassical model in Figure 5.6 predicts a continuous redshift of the LR for the off-resonant case (panels a), and increasing LR–UR splitting for the resonant case (panel b) with decreasing gap size D , consistent with the classical LRA predictions. As compared to these classical results, the main difference is that the semiclassical model shows a weakening and a broadening of the LR upon decreasing D , which points toward the role of nonlocal optical effects that can also affect the electromagnetic response for such a small system [ 336 ] as discussed in detail in Chapter 3. Therefore, the comparison of the results obtained within the classical LRA, the semiclassical model, and the TDDFT approach as employed in this chapter allows us to conclude that the blueshift of the LR for decreasing D below D≤ 26 a0 for the off-resonant QE–MNPs system, as well as the weakening of the LR–UR splitting for the resonant case, are only observed when the hybridization between the MNPs and the QE electronic orbitals is possible, i.e., when the (excited) electron can tunnel across the system. 5.2.3 Charge-transfer resonances at low frequencies Finally, we discuss in Figure 5.7 the role of the QE in triggering electron transport between the two MNPs across the junction in response to external illumination. In the last years, several works have identified the emergence of charge-transfer plasmons (CTP) supported by metallic vacuum junctions for gap separations typically below ∼ 0 . 4 nm and resonant frequencies of the order of a few electronvolts [ 43 , 68 , 83 , 84 , 86 , 342 , 357 ]. CTP are plasmonic resonances where a net electron transport occur between the MNPs that form the nanogap (see Subsection 2.3.2). 128 5.2. Results and discussion 12 a0 14 a0 16 a0 18 a0 isolated dimer 0=1 eV V V0=5 eV Figure 5.7:(a) Absorption spectra of the hybrid QE–MNPs system for low illumination frequencies ω = 0 − 1eV. Results are shown for a gap size D ranging from D = 12 a0 to D = 18 a0 , as indicated in the inset. The reference absorption spectrum of the isolated MNP dimer for D = 12 a0 is shown by the dashed black line. The situations for V0 = 1 eV (top) and V0 = 5 eV (bottom) in Eq. (5.3) are considered. (c) Color maps of the induced electron density (left) and the electron-current density along the z -direction (right) for an incident z -polarized electromagnetic plane wave of frequency ωCT = 0 . 11 eV. The gap distance is D = 16 a0 , and V0 = 1 eV . On the right-hand side panel, the boundaries of the jellium edges of the MNPs are indicated by dashed lines. The snapshots are taken at the instants of time when the absolute value of the total dipole moment (left-hand side panel) and of the electron-current density in the middle of the junction (right-hand side panel) are maximum. (c) Ground-state potential Veff along the symmetry z -axis of the hybrid QE–MNPs system for D= 16 a0. Moreover, as pointed out in previous works [ 107 – 109 , 329 , 353 ], the presence of a QE bridging a metallic nanogap substantially modifies the charge-transfer properties of the system and triggers out the emergence of low-frequency resonances associated with electron transport between the MNPs. In our study, the absorption cross-section σabs ( ω )shown in Figure 5.7a for two different QEs characterized by V0 = 1 eV (top) and V0 = 5 eV (bottom) reveals that a charge-transfer resonance emerges in the low-frequency region, ωCT ∼ 0 . 1 − 0 . 2 eV , for gap sizes D = 12 a0− 18 a0 ( D≈ 0 . 6 − 0 . 95 nm , thus larger than typical tunneling distances in metal–vacuum–metal junctions studied in Subsection 2.3.2). Our results are consistent with the findings reported in the literature [ 74 , 107 , 108 , 353 ]. This new resonance is only activated due to the presence of the QE (see the response of the isolated dimer depicted by the dashed line), and it blueshifts and strengthens considerably when decreasing interparticle 129 Chapter 5. Electronic exciton–plasmon coupling in a nanocavity beyond the electromagnetic interaction picture distance. The charge-transfer character of the mode is clearly revealed by the induced electron density shown in Figure 5.7b (left-hand side panel), with each MNP exhibiting a monopolar electron density pattern of opposite sign, and it is further corroborated by the electron-current density along the z -direction (righthand side panel), which clearly shows that electrons shuttle from one MNP to another. In the studied QE–MNPs system, the observed charge-transfer resonances at ωCT ∼ 0 . 1 − 0 . 2 eV emerge because for gap separations of D∼ 18 a0 and below, the QE gives rise to a decrease of the potential barrier close to the dimer axis below the Fermi level of the system (Figure 5.7c), so that even a classically-allowed over-the-barrier electron transport between the MNPs becomes possible. Thus, the low-frequency charge-transfer plasmon reported here can be understood as a consequence of the ballistic electron transport and does not require the tunneling mechanism aid by a localized state at the QE [109]. 5.3 Summary In summary, in this chapter we have identified the role played by electronic coupling in the optical response of a canonical hybrid system consisting in a twolevel quantum emitter (QE) placed in a nanogap formed by two spherical metal nanoparticles (MNPs). Using a fully TDDFT model for both the QE and the MNPs, we have demonstrated the quenching of the QE exciton originated by the hybridization of the excited states localized at the QE and the electronic states of the MNPs. This exciton quenching drastically affects the optoelectronic response of the hybrid QE–MNPs system for small gap separations. For example, it gives rise to a blueshift of the lower resonance (LR) with decreasing gap separation distance below D≲ 26 a0 for situations where the QE transition frequency ωQE is out of resonance with the main plasmonic mode of the MNP dimer, in contrast to the classical LRA calculations that predict a continuous redshift. Further, exciton quenching produced by the electronic interaction also leads to a drastic attenuation of the LR–UR splitting in resonant QE–MNPs systems, thus frustrating the strong coupling predicted by classical LRA simulations. On the other hand, depletion of the potential barrier within sub-nanometric gaps due to the presence of the QE gives rise to a low-frequency electron-transfer resonance at ωCT ∼ 0 . 2 eV , even for situations where the electronic states of the QE do not act as a gateway for electron transport between the MNPs. Our findings are expected to qualitatively apply for plasmon–exciton systems irrespective of the specific electronic structure of the nanoconstituents, since they are based on general and robust quantum-mechanical phenomena such as electron tunneling and electron transfer between the MNPs and the QE. Thus, the results obtained in this chapter stress the need to consider the QE–MNPs electronic coupling, in addition to the standard electromagnetic interaction, in order to unveil fundamental quantum effects related to charge transfer, often affecting practical implementation of nanoscale sources of photon emission and optoelectronic nanodevices. 130 Chapter 6 SECOND-HARMONIC GENERATION FROM A QUANTUM EMITTER COUPLED TO A METALLIC NANOPARTICLE In previous chapters, we focused on the analysis of quantum effects emerging in the optical and electronic response of plasmonic systems in situations where the intensity of the external illumination is weak and thus the response is linear. However, when the intensity of the external illumination is strong, the excitation of plasmonic resonances in metallic nanoparticles (MNPs) can also lead to nonlinear effects that can be useful for (bio-)imaging [ 358 – 360 ] or for generation of extreme-ultraviolet attosecond laser pulses [ 361 ], among others [ 36 , 37 ]. In particular, second-harmonic generation (SHG), whereby two photons at the fundamental frequency are absorbed to emit one photon at the second-harmonic frequency, is at the focus of very active research owing to its practical and fundamental interest [ 362 – 371 ]. In this context, it has been shown that plasmonic nanostructures resonant at the fundamental or at the second-harmonic frequency (or at both frequencies) can give rise to considerable enhancement of SHG [ 294 , 367 , 368 , 372 – 383 ]. Recent experiments have also shown the polarization-resolved probing of the nonlinear near-field distribution of metallic nanostructures by using doubly resonant plasmonic antennas [ 384 ]. To achieve SHG, however, the symmetry of the system needs to be considered. For example, we show in Section 1.3 and Section 2.4 that, for typical plane-wave incidence, SHG is forbidden from nanostructures that are centrosymmetric. This nonlinear response is thus very sensitive to the geometry of the system and to surface effects that may eventually break the symmetry constraints and lead to the emission of light at the second-harmonic frequency [37,68,154,385–388]. 131 Chapter 6. Second-harmonic generation from a quantum emitter coupled to a metallic nanoparticle In this chapter, we study SHG from a coupled system consisting of a quantum emitter (QE) placed in the vicinity of a spherical MNP [ 111 , 389 , 390 ], as scketched in Figure 6.1a. The small individual centrosymmetric MNP does not allow for second-harmonic emission into the far field, but it creates second-harmonic near fields in the proximity of the MNP surface. The presence of the QE lifts the symmetry constraints and allows for SHG. When the electronic transition frequency of the QE, ωQE , is resonant with the second harmonic of the incident frequency, the QE plays the role of an optical resonator, which efficiently couples to the nonlinear electric near field induced close to the MNP (see Section 2.4), transduces this near field into the far field, and thus produces SHG [ 384 ]. This QE–MNP system thus enables frequency conversion and allows for its control. To calculate the nonlinear response of the coupled system and to reveal the physical mechanisms behind SHG in this situation, we use TDDFT calculations [ 70 , 71 ] based on the wavepacket propagation (WPP) method introduced in Section 2.2. With the insights obtained from the TDDFT simulations, we develop a semi-analytical model that accurately reproduces the TDDFT results. This semi-analytical model also allows for addressing more general and complex situations beyond the reach of TDDFT, making possible a detailed study of the sensitivity of SHG to different parameters that characterize the system. In particular, we demonstrate the polarization conversion of the nonlinear signal, as well as the existence of various regimes of SHG determined by the intrinsic losses of the QE. The methodology and results presented in this chapter can pave the conceptual road for enhancing and optimizing SHG mediated by QEs coupled to plasmonic systems [391,392]. 6.1 System and methods We consider a QE located in the proximity of a spherical sodium MNP. The MNP is characterized as in Section 2.4 (Wigner–Seitz radius rs = 4 a0 , Ne = 1074 conduction electrons, and radius a = 40 . 96 a0 ). The dipolar plasmon (DP) resonance of the individual MNP is at ωDP = 3 . 17 eV , and the quadrupolar plasmon (QP) resonance at ωQP = 3 . 4 eV . A Gaussian-like external excitation Vext ( r, t ) given by Eq. (2.79) with fundamental frequency ω , duration σ = 5 × 2 π/ω , and intensity I0 = 10 10 W cm −2 (amplitude E0 = 4 . 8 × 10 −4au ) is used within TDDFT in this chapter. Importantly, as discussed in Section 2.4, such spherical MNP cannot emit second-harmonic light into the far field due to symmetry constraints, however second-harmonic fields with a quadrupolar pattern are induced in the proximity of the MNP because the inversion symmetry is locally broken at the surface (right-hand side panel in Figure 6.1b). The transition frequency of the QE is set to be resonant with the second harmonic of the fundamental frequency of the external illumination, ωQE = 2 ω , so that a variation of ω in our calculations implies simultaneous variation of ωQE . The QE plays the role of an optical resonator, sensitive to the second-harmonic electric near field [ 384 ]. We model the QE as a point-like dipole as described in Subsection 1.4.1, using an oscillator strength α0 QE = 1 au in Eq. (1.54) . The value of the intrinsic damping parameter γQE is 132 6.1. System and methods Figure 6.1:(a) Sketch of the system studied in this chapter: the radius of the spherical sodium MNP is a = 40 . 96 a0 ( ≈ 2 . 2nm), and the point-like QE is located at position rQE , at a distance d from the MNP surface. A Wigner–Seitz radius rs = 4 a0 is used to characterize the MNP within the jellium model. (b) Color maps of the real part of the radial component of the electric near field Eind ( r, ω )induced at the fundamental ( ω , left-hand side panel) and at the second-harmonic frequency (2 ω , right-hand side panel) by a z -polarized Gaussian electromagnetic pulse with fundamental frequency ω = 1 . 585 eV and intensity I0 = 10 10 W cm −2 incident at the individual spherical MNP in the absence of the QE (same results are also shown in Figure 2.8). Results are rotationally symmetric with respect to the z -axis, and they are shown in the ( x, z )-plane normalized to unity. Red and blue colors are used for positive and negative values, respectively (white for zero). varied in this chapter within the range γQE = 0.1eV −10−7eV. The expectation value of the QE dipole moment, pQE ( t ), evolves in time according to [22] [see Eq. (1.54)]: ¨ pQE(t) + γQE ˙ pQE(t) + ω2 QEpQE(t) = α0 QEEtot(rQE, t),(6.1) where the total electric field Etot ( rQE, t )acting on the QE position rQE is given by the sum of the incident laser pulse (with amplitude E0 , duration σ = 5 × 2 π/ω , and arrival time t0= 5σ), Eext(t) = ˆ zE0cos(ω(t−t0)) e−(t−t0 σ)2 ,(6.2) and the field Eind ( rQE, t )induced by the MNP [Eq. (2.34) ] at the position rQE of the QE, Etot ( rQE, t ) = Eext ( t ) + Eind ( rQE, t ). Note that Eind ( rQE, t )includes the reaction of the MNP not only to the incident pulse, but also to the electric field induced by the QE. It thus also accounts for the QE self-interaction (see Subsection 1.4.2). We describe in Appendix Chow we solve Eq. (6.1) in this thesis. The QE dipole pQE ( t )acts as a radiation source emitting into the far field as well as affecting the dynamics of the conduction electrons of the MNP. Because of the small size of the system, retardation effects can be neglected, so that the QE placed at a position rQE near the MNP creates an electrostatic potential given by [Eq. (3.3)]: VQE(r, t) = −pQE(t)·r−rQE |r−rQE|3.(6.3) 133 Chapter 6. Second-harmonic generation from a quantum emitter coupled to a metallic nanoparticle Thus, the Kohn–Sham Hamiltonian ˆ H [ n ( r, t )] within TDDFT [Eq. (2.19) ], and on the MNP electrons, is given by ˆ H[n(r, t)] = ˆ T+VH[n(r, t)] + Vxc[n(r, t)] + Vext(r, t) + VQE(r, t),(6.4) where, as discussed in Section 2.1 and Section 2.2, n ( r, t )is the time-dependent electron density [Eq. (2.20) ], ˆ T = −1 2∇2 is the kinetic-energy operator, VH [ n ( r, t )] is the Hartree potential [Eq. (2.23) ], Vxc [ n ( r, t )] is the exchange–correlation potential [Eq. (2.12) ] calculated using the kernel of Gunnarsson and Lundquist [ 198 ], and Vext ( r, t )is the external potential given by Eq. (2.79) that drives the QE–MNP system. We employ the WPP algorithm in spherical coordinates as described in Subsection 2.3.1. In this case, the time-dependent Kohn–Sham equations given by Eq. (2.19) and Eq. (2.20) are solved self-consistently together with Eq. (6.1) , Eq. (6.3) , and Eq. (6.4) . These equations are solved in time domain, and the time-to-frequency Fourier transform given by Eq. (2.80) is used to obtain the frequency-resolved quantities of interest such as the nonlinear dipole moment induced at the MNP pMNP (Ω), and at the QE pQE (Ω). The total dipole moment p (Ω) is given by the sum of both, p (Ω) = pMNP (Ω) + pQE (Ω). Here, the symbol Ωis used to denote the frequency of the induced dipole moments, since we are considering the nonlinear regime where ω = Ω in general. A Gaussian filter F ( t ) given by Eq. (2.81) is used in the Fourier transforms to partially account for decay and dephasing processes of the collective density oscillations that are not included in the present ALDA-TDDFT approach (see Section 2.2). 6.2 Results and discussion 6.2.1 TDDFT results We initially place the QE at the z -axis, corresponding to the direction of polarization of the incident laser pulse, at a distance d = 18 a0 ( ≈ 0 . 95 nm ) from the MNP surface. For this geometry, only z -polarized dipole moments are induced in the QE and in the MNP. The system then possesses cylindrical symmetry with respect to the z -axis, which greatly reduces the computational demands for the TDDFT calculations. The frequency of the external excitation ω = 1 . 585 eV is first considered such that its second harmonic matches the DP frequency of the MNP, 2 ω = ωDP = 3 . 17 eV , and the intrinsic damping parameter of the QE is set to γQE = 0.1eV. Figure 6.2 shows that the coupled QE–MNP system features strong emission at both odd and even harmonics, in contrast to the individual MNP that only emits at odd-harmonics due to the inversion symmetry of the system (see Section 2.4). In this figure, the intensity spectrum of the total induced dipole moment |p (Ω) |2 calculated within TDDFT for the coupled QE–MNP structure is shown by the blue line, revealing clear peaks at even harmonics Ω = 2 ω, 4 ω and 6 ω . The reference results obtained for the nonlinear response of the individual MNP (without QE) are shown by the dashed red line (only harmonics at Ω = ω, 3 ω, 5 ω and 7 ω are 134 6.2. Results and discussion Figure 6.2: Intensity spectrum of the total dipole moment |p (Ω) |2 of the coupled QE–MNP system (solid blue line) and of the individual MNP (dashed red line). Results are obtained for an incident z -polarized Gaussian electromagnetic pulse with fundamental frequency ω = 1 . 585 eV and intensity I0 = 10 10 W cm −2 . The QE is located at the z -axis, at a distance d = 18 a0 (≈0.95 nm) from the MNP surface, and it is characterized by a transition frequency ωQE = 2ω, an intrinsic damping parameter γQE = 0 . 1 eV , and oscillator strength α0 QE = 1 au , following Eq. (6.1) . Results are shown as a function of the frequency measured in units of the fundamental frequency ω . In the inset, the solid blue line corresponds to the same result as in the main figure (with QE), and the dashed black line corresponds to the results obtained for the transition frequency of the QE resonant with the fourth harmonic of the incident light, ωQE = 4 ω , with ω= 0.79 eV. observed, as discussed in Figure 2.7). The even harmonics in the far field from the coupled QE–MNP system emerge because the QE breaks the reflection symmetry with respect to the ( x, y )-plane, and thus the total inversion symmetry of the system [ 388 ]. Note that the spectra in Figure 6.2 are artificially broadened by application of the Gaussian filter given in Eq. (2.81) that allows for introducing losses in the system. We show in Figure 6.3 the effect of the Gaussian filter given by Eq. (2.81) . While the nonfiltered spectrum of the hybrid MNP-QE system (blue line) appears quite noisy, the filtered one presents very well-defined high-harmonic peaks (red line). All the peaks are broadened by the filter, which results in an attenuation of the maximum value of the peaks as compared to the nonfiltered signal. Figure 6.3 thus illustrates how this filtering procedure allows us to reach convergent spectral response at high-harmonic frequencies. In Subsection 6.2.2, we develop a semi-analytical method that allows us to overcome the difficulties of the ALDA-TDDFT calculations to incorporate losses [72,241,242]. The resonance between the transition frequency of the QE, ωQE , and the second harmonic of the incident pulse strongly enhances the intensity emitted by the system at 2 ω . To illustrate this resonance effect, we show in the inset of Figure 6.2 the results obtained for a different situation. The QE transition frequency ωQE in this case is set to be resonant with 4 ω (dashed black line), and the system is illuminated by a Gaussian pulse with fundamental frequency ω such that the fourth harmonic matches the frequency of the MNP dipolar plasmon, 4 ω = ωQE = ωDP = 3 . 17 eV. 135