<figure><img src="image_1.png" /></figure>

HAL
open science

# Computational Methods and Theoretical Models for the Study of Plasmonic and Hybrid Nanostructures

Nicolas Large

► To cite this version:

Nicolas Large. Computational Methods and Theoretical Models for the Study of Plasmonic and Hybrid Nanostructures. Materials Science [cond-mat.mtrl-sci]. Paul Sabatier Université Toulouse III - Paul Sabatier (UPS), 2023. tel-04181991

HAL Id: tel-04181991

https://laas.hal.science/tel-04181991

Submitted on 16 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Université Paul Sabatier – Toulouse III

<figure><img src="image_2.png" /></figure>

UNIVERSITÉ
TOULOUSE III
PAUL SABATIER

<figure><img src="image_3.png" /></figure>

Manuscrit

présenté pour l'obtention de

L'Habilitation à Diriger des Recherches

par

Nicolas Large

Department of Physics and Astronomy
The University of Texas at San Antonio

Computational Methods and Theoretical Models for
the Study of Plasmonic and Hybrid Nanostructures

Soutenance le 22 Juin 2023

devant le jury composé de:

**Rapporteurs:**
Renaud Bachelot
Jeremy Baumberg
Teri W. Odom
Professeur, Université de Technologies de Troyes
Professeur, University of Cambridge
Professeur, Northwestern University

**Examinateurs:**
Stéphane Calvez
Jennifer A. Dionne
Naomi J. Halas
HDR, Chargé de Recherche, LAAS-CNRS
Professeur, Stanford University
Professeur, Rice University

**Parrain:** Adnen Mlayah Professeur, Université Paul Sabatier – Toulouse III
# DIPLOME D'HABILITATION A DIRIGER DES RECHERCHES

## RAPPORT DE SOUTENANCE

Nom du candidat(e) : Dr. NICOLAS LARGE

Date de soutenance : 22 Juin 2023

Nom du Président du jury : Dr. RENAUD BACHELOT

Dr. Nicolas Large presented in English his research activities to an international jury through a clear and dynamic presentation, summarizing the highly coherent scientific work he has carried out since his Ph.D. in France and Spain until the development of his group, from 2016, in the department of Physics & Astronomy, University of Texas at San Antonio (UTSA). His presentation highlighted his scientific expertise in hybrid nanoplasmonics, as well as his ability to forge long-term collaborations, supervise young researchers, and obtain funding for numerous structuring projects. He also highlighted the richness of the subjects he tackled, which were conceptually complex and extremely innovative. They concern mainly new computational methods and original theoretical models but in strong connection with experimentalists. The rich discussions with the jury demonstrated Nicolas Large's vision and perspective on his field of research, as well as the interest of his scientific projects that are focused on high-potential subjects in many topical domains such as nanosensing, acousto-plasmonics, thermoplasmonics, metasurfaces, and magnetoplasmonics. He was able to discuss lucidly with the jury many of the problems and issues, and new directions that need to be developed in these areas.

The jury unanimously considered that Dr. Large's research works, responsibilities and supervision fully justify the award of the Habilitation à Diriger des Recherches of the Université of Toulouse, delivered by Paul Sabatier University.

## EMARGEMENT DES MEMBRES DU JURY

Dr. Renaud Bachelot

<figure><img src="image_4.png" /></figure>

Dr. Adnen Mlayah

<figure><img src="image_5.png" /></figure>

Dr. Jeremy Baumberg

<figure><img src="image_6.png" /></figure>

Dr. Teri W. Odom

<figure><img src="image_7.png" /></figure>

Dr. Stéphane Calvez

<figure><img src="image_8.png" /></figure>

Dr. Jennifer Dionne

<figure><img src="image_9.png" /></figure>

Dr. Naomi J. Halas

<figure><img src="image_10.png" /></figure>
*To my loving wife and wonderful children
who give me support and strength every day.*
<figure><img src="image_11.png" /></figure>
Acknowledgements

First and foremost I am extremely grateful to the amazing mentors I have had over the years, starting with Prof. Adnen Mlayah. Adnen has served as a mentor since my undergraduate studies (License) and his invaluable advice, continuous support, and patience throughout my career have forged the academic and scientist I am today. After 20 years, Adnen has become more than a mentor to me and his friendship remains a guiding force in my professional career. I would also like to thank Prof. Peter Nordlander who served as my first postdoc advisor. Peter has taught me all I know about Academia in the US. During the past 10 years he has shown me tremendous support and has helped me through challenging times. Peter remains a solid foundation of my academic career and one of the two persons I turn to when seeking advice. Adnen and Peter have shared their immense knowledge and experience countless times throughout the year, teaching me how to conduct cutting-edge research and publishing high-level research articles.

I would also like to express my deepest appreciation to the members of the committee for taking the time to read through the present manuscript and their invaluable support.

Thanks also to all my colleagues all around the world and with whom I collaborated over the years (in alphabetic order): Drs Carmen Afonso (Institute of Optics, Madrid), Daniel Aili (Linköping University), Javier Aizpurua (DIPC, CSIC), Denis Boudreau (Université Laval), Tapas Chini (Saha Institute of Nuclear Physics), Roland Coratger (CEMES, CNRS), Nien-Hui Ge (UC Irvine), José Gonzalo (Institute of Optics, Madrid), Naomi Halas (Rice University), Chen Hongyu, (NTU Singapore), Jean-Roch Huntzinger, (Université de Montpellier), Hao Jing (George Mason University), Bernard Jusserand (INSP), Todd Krauss (University of Rochester), Dmitry Kurouski (Texas A&M University), Stephan Link (Rice University), Jun Lou (Rice University), Kremena Makasheva (LAPLACE, CNRS), Katie Mayer (UTSA), Chad Mirkin (Northwestern University), Gesuri Morales-Luna (Universidad Iberoamericana), Otto Muskens (Universtiy of Southampton), Sina Najmaei (US Army Research Lab.), Junghwan Oh (Pukyong National University), Teri Odom (Northwestern University), Nancy Ornelas-Soto (Monterrey Tech.), Renaud Péchou (CEMES, CNRS), Arturo Ponce Pedraza (UTSA), Alejandro Reyes-Coronado (UNAM), Emilie Ringe (University of Cambridge), Lucien Saviot, Université de Bourgogne), George Schatz (Northwestern University), Jennifer Shumaker-Parry (University of Utah), Sudhiranjan Tripathy (IMRE, A\*STAR), Alex Urban (LMU Munich), Richard Van Duyne† (Northwestern University), Rémi Vincent, (UTT), Hui Wang (University of South Carolina), and Shan Wang (Stanford University).

Finally, I would like to express my gratitude to my wife, Juan, and my children, Logan and Alison. They have born with patience and understanding my late-night work, long working weekends, and frequent travels. Their Love, smiles, and tremendous support strongly contributed to my professional success.
Abstract

This manuscript synthesizes my scientific activity since Ph.D.; covering the period from
November 2011 to June 2023. While some parts of the work presented here were first carrier
out in the context of my two postdoctoral positions at Rice University (2012-2014) and
Northwestern University (2014-2016), the bulk of this scholarly work has been performed
during my time as a faculty in the Department of Physics and Astronomy at the University of
Texas at San Antonio (UTSA), where I am currently an Associate Professor with Tenure.

For this manuscript, I have chosen to present selected pieces of my research activity on the
development of novel theoretical approaches and computational methods for the study of the
optical, thermo-optical, and photothermal properties of nanomaterials. A general synthesis of
my scientific activity is provided in Chapter 1. It gives a global overview of my career
trajectory, research projects and collaborations, and student advising activity since I joined
UTSA. Chapters 2-5 then discuss in greater details several significant research projects which
include: (i) a new electron energy-loss spectroscopy simulation method, (ii) a new method for
the calculation of photonic band structures in finite periodic systems, (iii) a new plasmon-
vibration interaction theoretical framework, and (iv) a new plasmonic heating simulation
method. A final chapter will provide some scientific perspectives and general research
directions for the near, mid, and long terms. Lastly, I will wrap-up with a summary of my
contribution in *computational methods and theoretical models for the study of plasmonic and
hybrid nanostructures*.

**KEYWORDS:** Nanophotonics; Plasmonics; Condensed Matter Theory; Computational Physics; Electrodynamics; Nanomaterials; Optical Materials.
# Table of Contents

<table><tr><td>Acknowledgements</td><td>iii</td></tr><tr><td>Abstract</td><td>iv</td></tr><tr><td>Table of Contents</td><td>v</td></tr><tr><td>List of Figures</td><td>vii</td></tr><tr><td>List of Tables</td><td>x</td></tr><tr><td>Introduction</td><td>1</td></tr><tr><td>1 | General Synthesis</td><td>3</td></tr><tr><td>A. Career Overview</td><td>3</td></tr><tr><td>B. Research Projects Overview</td><td>4</td></tr><tr><td>1. Elementary Excitation Interactions</td><td>4</td></tr><tr><td>2. Optical, Photothermal, and Thermo-Optical Properties of Nanomaterials</td><td>5</td></tr><tr><td>3. Plasmonics for Chemistry and Biomedicine</td><td>5</td></tr><tr><td>C. Student Supervision and Mentoring</td><td>6</td></tr><tr><td>D. Funded Research Projects</td><td>6</td></tr><tr><td>2 | Electron Energy-Loss Spectroscopy in Finite-Difference Time-Domain</td><td>9</td></tr><tr><td>A. Related Publications</td><td>9</td></tr><tr><td>B. Background</td><td>9</td></tr><tr><td>C. Method</td><td>10</td></tr><tr><td>D. Results and Discussion</td><td>12</td></tr><tr><td>1. Isolated Nanostructure</td><td>12</td></tr><tr><td>2. Interacting Nano-Objects</td><td>15</td></tr><tr><td>3. Supported Nanostructures</td><td>16</td></tr><tr><td>4. Complex Nanostructures</td><td>18</td></tr><tr><td>5. Periodic Nanoparticle Arrays</td><td>21</td></tr><tr><td>E. Conclusions</td><td>22</td></tr><tr><td>3 | Photonic Band Structures of Finite Periodic Nanostructures</td><td>23</td></tr><tr><td>A. Related Publications</td><td>23</td></tr><tr><td>B. Background</td><td>23</td></tr><tr><td>C. Methods</td><td>25</td></tr><tr><td>1. Infinite Case</td><td>25</td></tr><tr><td>2. Finite Case</td><td>26</td></tr><tr><td>3. Effective Medium Model</td><td>27</td></tr><tr><td>D. Results and Discussion</td><td>27</td></tr><tr><td>E. Conclusions</td><td>33</td></tr></table>
<table><tr><td>4</td><td>Raman Energy Density in the Context of Acousto-plasmonics</td><td>35</td></tr><tr><td>A.</td><td>Related Publications</td><td>35</td></tr><tr><td>B.</td><td>Background</td><td>35</td></tr><tr><td>C.</td><td>Theoretical Framework</td><td>36</td></tr><tr><td>1.</td><td>Raman Energy Density (RED)</td><td>36</td></tr><tr><td>2.</td><td>Raman Scattering Spectrum</td><td>38</td></tr><tr><td>D.</td><td>Model and Methods</td><td>39</td></tr><tr><td>1.</td><td>Computational Model</td><td>39</td></tr><tr><td>2.</td><td>Elastodynamics: Resonant Ultrasound Spectroscopy (RUS) Method</td><td>39</td></tr><tr><td>3.</td><td>Electrodynamics: Discrete Dipole Approximation (DDA)</td><td>39</td></tr><tr><td>E.</td><td>Results and Discussion</td><td>40</td></tr><tr><td>1.</td><td>Acoustic Vibrations</td><td>40</td></tr><tr><td>2.</td><td>Plasmonic Near-Field</td><td>42</td></tr><tr><td>3.</td><td>Raman Energy Density (RED)</td><td>43</td></tr><tr><td>4.</td><td>Acoustic Raman Spectra</td><td>45</td></tr><tr><td>5.</td><td>Raman Selection Rules</td><td>46</td></tr><tr><td>F.</td><td>Conclusions</td><td>48</td></tr><tr><td>5</td><td>Plasmonic Photothermal Heating Effects in Nanostructures</td><td>49</td></tr><tr><td>A.</td><td>Related Publications</td><td>49</td></tr><tr><td>B.</td><td>Background</td><td>49</td></tr><tr><td>C.</td><td>Method</td><td>50</td></tr><tr><td>1.</td><td>Optical Absorbance</td><td>50</td></tr><tr><td>2.</td><td>Absorbed Power Density</td><td>51</td></tr><tr><td>3.</td><td>Heat Transfer Calculation</td><td>52</td></tr><tr><td>D.</td><td>Results and Discussion</td><td>52</td></tr><tr><td>1.</td><td>Single Au Nanoparticles</td><td>52</td></tr><tr><td>2.</td><td>Periodic Au Nanoparticle Arrangements</td><td>55</td></tr><tr><td>3.</td><td>Random Ag Nanoparticle Ensembles</td><td>57</td></tr><tr><td>4.</td><td>TERS Tip</td><td>58</td></tr><tr><td>E.</td><td>Conclusions</td><td>59</td></tr><tr><td>6</td><td>Perspectives</td><td>61</td></tr><tr><td>A.</td><td>Acousto-Plasmonics</td><td>61</td></tr><tr><td>1.</td><td>Application of the RED Framework to other Optomechanical effects</td><td>61</td></tr><tr><td>2.</td><td>Optical Forces Effect on Acoustic Vibrations</td><td>62</td></tr><tr><td>B.</td><td>Multidimensional Periodic Systems</td><td>63</td></tr><tr><td>1.</td><td>Surface Lattice Resonances in Plasmonic Periodic Structures</td><td>63</td></tr><tr><td>2.</td><td>New Polaritonic Modes in Hybrid Plectonic Periodic Structures</td><td>64</td></tr><tr><td>C.</td><td>Photothermal and Thermo-Optical Properties</td><td>65</td></tr><tr><td>1.</td><td>Temperature Effects on the Plasmon-Exciton Coupling</td><td>65</td></tr><tr><td>2.</td><td>Temperature Effects on the Plasmon-Phonon Coupling</td><td>66</td></tr><tr><td>Summary</td><td></td><td>67</td></tr><tr><td>CV</td><td></td><td>69</td></tr><tr><td>References</td><td></td><td>77</td></tr></table>
**Figure 7.** (a) STEM-HAADF images of a hollow AuAg nanorod and experimental EELS maps of each of the four LSPRs. (b) Decomposition of the calculated loss probability [$\Gamma_{EELS}$, Eq. (8)] maps calculated using the EELS-FDTD method for the center state nanorod (i.e., $\Delta = 0$ nm). (c) Calculated surface charge density distribution (red: positive charge, blue: negative charge) for each of the four LSP modes obtained in the FDTD calculations from panel b. 21

**Figure 8.** Left: Loss probability spectra calculated for a 6×6 NP array. The EELS spectra are color-coded with the various impact parameters shown in the inset. The extinction cross-section calculated for a plane wave excitation is shown as a comparison (blue). Right: Loss probability maps and surface charge density distributions calculated at 4 different energies. 22

**Figure 9.** Approaches to calculate the properties of photonic crystals, i.e., infinite lattice for band structure calculations, effective medium theory for optical properties like reflectance, and proposed FDTD simulation for finite supercrystals. 24

**Figure 10.** (a) Real space simple cubic lattice (left) formed with Au NPs showing the dipole sources (blue arrows) and field monitors (red crosses) within the unit cell used for the infinite case simulation and the first Brillouin zone of the reciprocal cubic lattice (right) showing the high symmetric path ($\Gamma$-X-M-R-Γ) used for the band structure calculations. (b) Block condition. (c) Supercrystals with cubic, rhombic dodecahedral, and spherical habits, used in the finite case simulations made up of Au NPs. (d) Simulation layout of the finite supercrystal (9 unit cells are shown) consisting in Au NPs (yellow spheres), electric dipoles (blue arrows), and field monitors (red crosses). 27

**Figure 11.** Photonic band structures of an infinite cubic lattice composed of 80 nm Au NPs and with a lattice parameter of (a) 100 nm and (b) 130 nm. 28

**Figure 12.** Photonic band structures of the finite cubic-lattice supercrystals with (a,b) cubic, (c,d) rhombic dodecahedral, and (e,f) spherical habits. The supercrystals are composed of 80 nm Au NPs with lattice parameters of 100 nm (left panels) and 130 nm (right panels). 29

**Figure 13.** Mode splitting to plasmon mode energy ratios ($\Omega R/\omega 0$) at the X point for the finite supercrystals with (a,b) cubic and (c,d) rhombic dodecahedral habits $\hbar\omega 0 \sim 2.3$ eV. 30

**Figure 14.** Reflectance at the X point ($k = kX$, left panels) with the FPMs indicated by red markers (the red arrow in the inset indicates the incident wavevector). Dispersion relation (right panels) of Au NPs supercrystals with a lattice parameter of 100 nm. (a,b) Cubic, (c,d) rhombic dodecahedral, and (e,f) spherical habits. In the dispersion relations are included the light dispersion in water (cyan solid line), dispersion relation calculated with FDTD of a slab with effective refractive index of the supercrystal (white dashed line), and dispersion relation of the light in the effective medium (blue dashed line) with the FPMs obtained from the reflectance spectrum and applying the resonant condition of a Fabry-Pérot cavity, to find the k value, $kN = \omega nReNeff\omega/c$ (magenta crosses). 31

**Figure 15.** (a) Reflectance spectra calculated for an incident polarization along the nanorods long axis (blue) and short axis (red) of a finite supercrystal with cubic habit composed of Au nanorods with a length of 80 nm and a width of 40 nm, arranged in a simple cubic lattice with a lattice parameter of 100 nm. (b) Photonic band structure of the same finite cubic supercrystal. 33

**Figure 16.** Calculated absorption spectrum of a 5 nm AuNP in water showing the dipolar LSPR at 520 nm. Inset: 3D color map and (001) cross-sectional plane of the ENF at the LSPR. 40

**Figure 17.** (a) AuNP at rest. Top: NP deformed by the isotropic breathing $S001$ ($d_0$) mode (b) and the fivefold degenerated quadrupole $S2m1$ ($d_1 - d_5$) modes (c-g). Bottom: NP deformed by the anisotropic breathing $A_{1g}$ mode (h), the quadrupole $E_g$ [$E_{g(a)}$ and $E_{g(b)}$] modes (i-j), and the quadrupole $T_{2g}$ [$T_{2g(a)}$, $T_{2g(b)}$, and $T_{2g(c)}$] modes (k-m). 41

**Figure 18.** ENF induced by the dipolar LSP and modulated by the isotropic (top) and anisotropic (bottom) vibration modes. The cross-sections are taken where the maximum deformation is located; the (xyz) planes are indicated in each panel. (a) NP at rest, (b) isotropic breathing $S001$ ($d_0$) mode, (c-g) fivefold degenerated isotropic quadrupoles $S201$ ($d_1 - d_5$), (h) anisotropic breathing $A_{1g}$ mode, (i,j) anisotropic quadrupole $E_g$ [$E_{g(a)}$ and $E_{g(b)}$] modes, and (k-m) anisotropic quadrupole $T_{2g}$ [$T_{2g(a)}$, $T_{2g(b)}$],
Introduction

Multifunctional nanomaterials are of high interest to a wide range of technological applications ranging from optoelectronics to biomedicine, energy, sustainability, photonic devices, and applied chemistry. In order to develop multifunctional nanomaterials capable of addressing specific functions for targeted applications, it is crucial to first understand the fundamental physical properties and mechanisms involved in such complex nanostructures. Well-known examples are excitons, sustained in low-dimensional semiconductors, and localized surface plasmon resonances (LSPRs), sustained by nanostructures composed of noble metals, which are two types of collective electronic excitations whose spectra are in the visible range. Since their theoretical predictions by Rufus Ritchie in 1957,¹ localized surface plasmons continue to serve as robust foundations for nanotechnologies. Their properties are well known and their control has been used in various applications in photochemistry,²⁻⁵ enhanced spectroscopies,⁴,⁶⁻¹⁰ medicine,¹¹ ¹² imaging,¹³ sensing,¹⁴⁻¹⁷ active light manipulation/circuitry,¹⁸⁻²² and solar applications.²⁰, ²³, ²⁴

In the past decade, a great deal of effort has been made to improve the performance of plasmonic-based applications. Adding multiple levels of functionality to a nanoparticle (NP) can be done in various ways, such as doping or nanomaterial structuration. The latter can be obtained using multi-layered core-shell NPs where the core and shells perform different functions.¹⁹, ²⁵ Another route to tailoring novel optical properties in complex plasmonic-based structures consists of structuring the building blocks in a well-organized manner, forming one-dimensional (1D) chains,²⁶⁻²⁸ two-dimensional (2D) arrays,²⁸, ²⁹ or three-dimensional (3D) superlattices.²⁸, ³⁰, ³¹ Although there have been a multitude of studies on the influence of dimensionality and arrangements of polarizable particles, almost all of these studies have focused exclusively on organized monomaterial hyperstructures (e.g., pure plasmonic, pure dielectric).²⁸, ²⁹, ³² Hybrid multidimensional arrangements composed of multiple materials remain vastly unexplored with the exception of some recent works on 2D plasmonic transition metal dichalcogenides (TMDs) nanodisk arrays³³, ³⁴ and on 3D magneto-optical superlattices.³⁵ Another fundamental aspect resulting from interfacing dissimilar materials and nanostructuring is the emergence of interactions between LSPs and other type of elementary excitations including polaritons,³⁶ magnons,³⁷ and excitons.³⁸ Control of these interactions has also been used as a strategy to improve the functionality and performance of photonic applications. More importantly, each new elementary-excitation interaction has led to the emergence of new fundamental concepts such as plasmon hybridization, plexcitons, lattice plasmons, and ferroplasmons.³⁹

This manuscript contains a synthesis of my research on the computational and theoretical study of the optical properties of nanomaterials. As I will discuss, this encompasses the development of novel theoretical approaches and computational methods for the modeling of the properties of optical nanomaterials composed of semiconductors, plasmonic and catalytic metals, biomaterials, dielectric materials, and their hybrids. Examples of such developments,
which will be discussed in this manuscript, include new models and methods for the description of acoustic-vibration Raman scattering, electron-based spectroscopies, plasmonic heating, and photonic band structure calculations. These new numerical methods and theoretical models are then used to study nanostructured systems and in particular the effects and properties arising from (i) interactions between elementary excitations in hybrid and periodic nanostructures, (ii) finite-size effects in multi-dimensional periodic arrangements, and (iii) photothermal processes.

First, I will present a general synthesis of my scientific activity in Chapter 1. This will serve as a preamble to the other technical chapters and will place my work in a broader context. The first section provides a general overview of my academic carrier, starting with my Ph.D. work. The second section then discusses my three main research trusts and describes the different scientific projects and collaborations in which I have been involved as an independent researcher. Note that, due to the theoretical nature of my research expertise, a significant part of my research activity is carried out in collaboration with experimental groups. Finally, the third and final section provides a short presentation of the students I have formally supervised in my research group between 2016 and 2023, and who have done all the heavy lifting on the research projects.

Second, Chapters 2-5 present in great details four significant projects illustrative of my contribution to the field of nanophotonics. Each of these chapters follows the same structures. First, I provide a list of related articles, where the students I supervised are underlined in the author list. Second, a short background section will provide some context for the work discussed in the chapter and is followed by a technical section describing the method we implemented. Finally, I discuss in details the results obtained with this method and summarize with some conclusions. A final chapter will provide some scientific perspectives and general research directions for the near, mid, and long terms in the main areas discussed in the 4 previous chapters.

Finally, I present a summary of my research activity, career trajectory, and contribution to the field of nanophotonics. This final summary aims to convince the reader that I have developed a successful, coherent, and well-thought research strategy by fostering research collaborations, supervising and mentoring students, and securing external extramural funding to support my research program.
1 |

General Synthesis

**A. Career Overview**

After obtaining a *License de Physique* (2005) and *Master Physique de la Matière* (2007) at the Université Paul Sabatier – Toulouse III, I obtained, in 2011, a dual Ph.D. degree in *Nanophysique* from the Université de Toulouse and in *Physics of Nanostructures and Advanced Materials* from the University of Basque Country in San Sebastián, Spain. During this time, I jointly worked at the Center of Materials Physics (CFM-CSIC) and Donostia International Physics Center (DIPC) in San Sebastián, and at the Centre d'Élaboration de Matériaux et d'Etudes Structurales (CEMES-CNRS) in Toulouse. Under the supervision of Drs Adnen Mlayah and Javier Aizpurua, I developed several novel approaches and numerical methods for the description of acoustic-phonon Raman scattering in semiconducting and metallic low-dimensional systems. My doctoral work led to the publication of 6 peer-reviewed articles and 2 conference proceedings.

From 2012 to 2014, I worked as a postdoctoral researcher in the group of Dr. Peter Nordlander, in the Department of Electrical and Computer Engineering at Rice University, Houston. There, I further developed my research expertise in the area of computational plasmonics, collaborating with experimental research groups both locally, nationally, and internationally as discussed in the following section. Later, from 2014 to 2016, I served as a Senior Research Fellow in the group of Dr. George Schatz, in the Department of Chemistry at Northwestern University, Chicago. During this second postdoc, I applied my computational and theoretical expertise to Chemistry problems and in particular to the modeling of optical spectroscopy techniques such as surface-enhanced and tip-enhanced Raman spectroscopies as well as second harmonic generation. My postdoctoral work led to the publication of 16 peer-reviewed articles.

Finally, in 2016, I was recruited as a tenure-track Assistant Professor in the Department of Physics and Astronomy at the University of Texas at San Antonio (UTSA). My recruitment was part of a program to strengthen the materials science program UTSA. The institution has an initiative that focuses on growing research in various areas, to include Advanced Materials and Materials Chemistry. Since 2016, I have built my own independent research group with extramural research funding (totaling $1.3M) and developed my own research program, which focuses on the computational and theoretical study of the optical properties of nanomaterials composed of semiconductors, plasmonic and catalytic metals, biomaterials, dielectric materials, and their hybrids. In 2022, as a result of my contribution and performance in Research, Teaching, and Service, I was promoted to Associate Professor with Tenure. Later, in the spring 2022, I joined the CEMES as a NanoX Visiting Researcher, where I conducted collaborative research for 3 months. To date, my independent work as a faculty has resulted in 20 peer-reviewed articles, 3 conference proceedings, and 22 oral (including 16 invited) presentations at international conferences.
## B. Research Projects Overview

As discussed above, my research program is centered on the computational and theoretical study of the optical properties of nanomaterials. This encompasses the development of novel theoretical approaches and computational methods for the modeling of the properties of optical nanomaterials composed of semiconductors, plasmonic and catalytic metals, biomaterials, dielectric materials and their hybrids. Examples of such developments include new models and methods for the description of acoustic vibration Raman scattering, electron-based spectroscopies, plasmonic heating, and photonic band structure calculations. My group also conducts a wide range of fundamental and applied studies in close collaboration with experimental groups. Although my research projects are intertwined, they are organized in three major thrusts:

### 1. Elementary Excitation Interactions

Hybrid nanomaterials, which are composed of dissimilar materials, sustain a variety of elementary excitations (also known as quasiparticles) such as excitons, plasmons, phonons, and polaritons, which interact. Understanding these interactions and developing new interaction mechanisms is crucial to understand and control the hybrid nanomaterials' physical properties, ultimately leading to the development of high accuracy quantum sensors, light emitting devices, new display technologies, and high-density data storage devices.

**Plasmonic-Excitonic Nanomaterials:** My group has contributed to the fundamental understanding of the plasmon-exciton interaction in systems composed of noble metal and semiconducting materials such as transition metal dichalcogenides (TMD)$^{33,40}$ and molecular J-aggregates.$^{38}$ In particular, we have reported a unique contribution of the TMD's exciton to plasmon-enhanced optical spectroscopies such as surface-enhanced Raman scattering (SERS)$^{33}$ and photoluminescence (PL).$^{40}$ These projects, done in collaboration with Drs. Péchou, Coratger, and Mlayah (CEMES), are the first combined theoretical experimental studies of plasmon-exciton mediated SERS and PL. This project has been supported by two Chateaubriand Fellowships from the Embassy of France in the U.S. ($24k in total) and a federal grant from the Office of Naval Research of the U.S. Department of Navy ($40k).

**Plasmonic-Photonic Periodic Nanomaterials:** My group has been investigating the optical properties of periodic arrangements of nanomaterials, which exhibit very novel properties as compared to their non-periodic counterparts. Such systems include 2D arrays of silicon pillars, gold/silver/aluminum nanoparticles, and nanoholes, as well as 3D gold nanoparticle supercrystals. Several projects involved collaborations with experimental groups, including Drs. Vincent (UTT), Chiang (U. Houston), and Mayer (UTSA). The most significant contribution from my group is the development of a novel computational method for the computation of photonic band structures in finite periodic arrays.$^{41}$ Our method, which is the first implementation for finite-size systems allows for unravelling effects occurring as a result of the finite size such as photonic cavity modes (Fabry-Pérot, whispering gallery), surface and edge effects, crystal habit effects, and changes in photonic mode lifetime, which all occur in realistic, experimentally fabricated systems. I will discuss this particular project in details in Chapter 3. Part of this project has received support through the Thomas Jefferson Fund from the FACE Foundation ($20k) and from the National Science Foundation ($60k).

**Plasmonic-Phononic Nanosystems:** My group theoretically investigates how acoustic vibrations (phonons) in nanomaterials modulate the surface plasmons in these nanostructures and therefore, how these vibrations influence their optical properties. We have developed a new mathematical formalism and computational method for the modeling of these effects in the framework of acoustic-phonon Raman scattering. This has recently led to the introduction of a new physical quantity, the Raman energy density (RED), which provides new insight
into these effects. Our newly developed method is the first non-atomistic method that allows for computing acoustic Raman spectra. I will discuss this particular project in details in Chapter 4. This research project has received major federal funding from the Office of Naval Research of the U.S. Department of Navy ($435k).

## 2. Optical, Photothermal, and Thermo-Optical Properties of Nanomaterials

When nanomaterials are optically excited, various physical phenomena can occur. In semiconductors, this excitation can lead to the generation of free charges in the conduction band, which leads to drastic changes of their optical properties. In metals, part of the electromagnetic energy is converted into heat, which also impacts their optical properties. Developing new models and methods to rigorously describe these light-matter interaction processes is crucial for the understanding and control of the nanomaterials optical properties.

Over the past few years, my group has been investigating the optical, photothermal, and thermo-optical properties of various nanomaterials. The major contribution in this research area has been the collaboration with the group of Dr. Mayer for the development of a complete computational approach for the rigorous modeling of photothermal processes.⁴² I will discuss this particular activity in details in Chapter 5. This research project has received major federal funding from the Army Research Office of the U.S. Department of Defense ($650k).

We have also shown that plasmon induced temperature increase has significant impact on the nanomaterial surrounding environment. In particular, we showed that (i) it can be used to trigger the contraction of thermo-responsive polymers (with Dr. Ye, UTSA),⁴³ (ii) it is responsible for molecular degradation in TERS experiments (with Dr. Kurouski, TAMU),⁴⁴ and (iii) it can be used in hyperthermia for cancer therapy (with Dr. Reyes-Coronado, UNAM and Dr. Oh, PNU).¹¹, ⁴⁵

However, the computational work carried out in my group has also shown that photothermal processes do not always contribute to plasmon-driven catalysis and that other processes such as hot electron generation are at play.⁴⁶, ⁴⁷

Finally, we recently showed that hybrid gold-silicon-gold nanodisk stacks can sustain higher-order localized surface plasmon modes at the gold-silicon interfaces.⁴⁸ These modes can be dynamically manipulated by changing the concentration of optically-excited free carriers in the silicon layer. This opens tremendous opportunities for the design of active photonic devices such as ultrafast optical switches.

## 3. Plasmonics for Chemistry and Biomedicine

When plasmons are optically excited in metallic nanostructures, a cascade of physical phenomena can be used to perform chemistry or interact with the environment. These include the generation of (i) strong localized electric fields, (ii) hot electrons, and (iii) heat. Understanding and controlling these effects allows for improving the efficiency of catalytic reactions important to our society, detecting environmental contaminants, as well as for drug delivery and theranostics.

Through various collaborations with experimental groups, I have actively contributed to the field of plasmon-mediated chemistry and biomedicine. We have provided significant insight into the contributions of photothermal effects and hot electrons in various catalytic reactions. This was rendered possible thanks to the development of a self-consistent computational method for the rigorous modeling of photothermal effects in single NPs and NP ensembles (cf. Section 2). My group has collaborated with Dr. Kurouski (TAMU) to show that photothermal effects were not responsible for the enhancement in catalytic activity observed in certain chemical reactions. Very recently, the accurate modeling of plasmonic electric field
done in my group has unraveled the critical role of rectified electric fields and Stark effect in the catalytic activity of bimetallic nanostructures.⁴⁶ In addition, my group has routinely collaborated with various experimental groups, such as Drs. Makasheva (LAPLACE-CNRS),⁴⁹ Mlayah (CEMES-CNRS),³³ Kurouski (TAMU),⁸,⁵⁰ and Ornelas-Sotto (TEC)⁵¹ to design and optimize SERS substrates and interpret SERS/TERS experiments. Part of this research effort has received support from the ConTex program from Mexico’s National Council of Science and Technology (CONACyT).

### C. Student Supervision and Mentoring

Since my hiring as an Assistant Professor in 2016, I have supervised and mentored a total of 37 students, including 21 undergraduate (UG)ª students from Physics, Computer Science, Mathematics, and Chemistry majors, 4 visiting students (1 B.S., 3 M.S., and 1 Ph.D. from France and Mexico), 10 graduate (GR)ᵇ students, and 1 postdoctoral researcher.

All of these students have been assigned at least one research project, have contributed to the science presented in this manuscript, and have been financially supported through various research assistantships, fellowships, and research grants. More details on my Advising and Mentoring activity can be found in the CV section on page 75. My students are listed as co-authors on 18 of the 23 articles I have published between 2016 and 2023. Additionally, my postdoc has published

one peer-reviewed article while under my supervision (J. Phys. Chem. C 2019, 123, 18), and for which I chose not to be listed as a co-author to demonstrate his independence in research.

As of June 2023, 5 graduate students have successfully defended their Ph.D. thesis [G. Morales Luna (2019), E. Peñafiel (2021), V. Nooshnad (2021), J. Rigor (2022), P. Lopez (2023)] and 2 undergraduate students have successfully defended their Honors thesisᶜ [J.I. Garcia (2019), R. Kraja (2022)]. Out of the 5 students who obtained their Ph.D. under my supervision, 4 of them have secured a permanent job (Morales Luna, Peñafiel, Rigor, Lopez) and 1 has started a second Ph.D. degree (Nooshnab).

<figure><img src="image_12.png" /><figcaption>Dr. Large's research group (2019)</figcaption></figure>

### D. Funded Research Projects

During the past 7 years, I had 12 projects which have been funded by various U.S. agencies, private foundations, and other sponsors, for a total of $1.3M. Out of these 12 projects, I served as the Principal Investigator (PI; i.e., project leader) on 7 projects and as co-PI on the other 5 projects. Eight of these projects have been successfully completed and have led to publications and conference presentations; the other 4 projects are still active to date.

Below, I present a summary of these funded projects, which include the list of co-PIs on the project with relative contributions, project title, funding agency, funding amount, and project period.

<footnote>ª UG is 4 years, meaning that Years 1-3 correspond to the *License*, while Year 4 corresponds to *Master 1*.</footnote>

<footnote>ᵇ All GR students are doctoral (Ph.D.) students.</footnote>

<footnote>ᶜ Equivalent to France's Master thesis</footnote>
**Project Title:** Acousto-plasmonic interaction for unconventional nanoscale light manipulation and integrated high accuracy sensors.

**PIs:** N. Large (PI, 100%)

**Sponsor:** Office of Naval Research (ONR), U.S. Department of the Navy (DoN)

**Funded Amount:** $434,446

**Project Period:** 7/2021 – 6/2024

**Project Title:** NSF Center for quantum electrodynamics for selective transformations.

**PIs:** N. Large (PI, 60%), Z.J. Tonzetich (co-PI, 40%)

**Sponsor:** National Science Foundation (NSF)

**Funded Amount:** $59,997

**Project Period:** 9/2021 – 8/2024)

**Project Title:** Nano-colloids of core-shell particles as a new approach for the optical sensing of contaminants of emerging concern.

**PIs:** K. Mayer (PI, 34%), N. Large (co-PI, 33%), N. Ornelas-Soto (co-PI, 33%)

**Sponsor:** ConTex Consortium

**Funded Amount:** $70,000

**Project Period:** 9/2020 – 8/2023

**Project Title:** HYMNE: unravel novel polaritonic states in hybrid metasurfaces and nanomaterials.

**PIs:** N. Large (PI, 50%), R.P. Vincent (co-PI, 50%)

**Sponsor:** FACE Foundation / Thomas Jefferson Fund

**Funded Amount:** $20,000

**Project Period:** 9/2019 – 8/2023

**Project Title:** Computational modeling of the thermo-optical and photothermal properties of hybrid plasmonic-excitonic nanostructures.

**PIs:** N. Large (PI, 100%)

**Sponsor:** NanoX

**Funded Amount:** €8,000

**Project Period:** 3/2022 – 4/2022

**Project Title:** Theoretical investigation of the optical and photothermal properties of plasmonic-based nanostructures.

**PIs:** N. Large (PI, 50%), A. Mlayah (co-PI, 50%)

**Sponsor:** Embassy of France in the U.S.

**Funded Amount:** $12,000

**Project Period:** 9/2019 – 8/2020
**Project Title:** Optical properties of plexcitonic nanostructures.

**PIs:** N. Large (PI, 50%), A. Mlayah (co-PI, 50%)

**Sponsor:** Embassy of France in the U.S.

**Funded Amount:** $12,000

**Project Period:** 9/2020 – 8/2021

**Project Title:** Hybrid plasmonics at UTSA: Investigating plasmonic/magnetic and plasmonic/biomolecular systems.

**PIs:** K. Mayer (PI, 30%), N. Large (co-PI, 25%), A. Ponce Pedraza (co-PI, 15%), K.L. Nash (co-PI, 15%), L. Brancaleon (co-PI, 15%)

**Sponsor:** Army Research Office (ARO), U.S. Department of Defense (DoD)

**Funded Amount:** $599,989

**Project Period:** 9/2018 – 9/2021

**Project Title:** Theoretical modeling of magneto-plasmonic effects in electrodynamic calculations.

**PIs:** N. Large (PI, 70%), A. Reyes Coronado (co-PI, 30%)

**Sponsor:** American Physical Society (APS)

**Funded Amount:** $2,000

**Project Period:** 3/2019 – 2/2024

**Project Title:** HSAP/URAP supplement: Hybrid plasmonics at UTSA.

**PIs:** K. Mayer (PI, 20%), N. Large (co-PI, 20%), A. Ponce Pedraza (co-PI, 20%), K.L. Nash (co-PI, 20%), L. Brancaleon (co-PI, 20%)

**Sponsor:** Army Research Office (ARO), U.S. Department of Defense (DoD)

**Funded Amount:** $22,500

**Project Period:** 6/2021 – 8/2021

**Project Title:** HSAP/URAP supplement: Hybrid plasmonics at UTSA.

**PIs:** K. Mayer (PI, 20%), N. Large (co-PI, 20%), A. Ponce Pedraza (co-PI, 20%), K.L. Nash (co-PI, 20%), L. Brancaleon (co-PI, 20%)

**Sponsor:** Army Research Office (ARO), U.S. Department of Defense (DoD)

**Funded Amount:** $25,000

**Project Period:** 6/2020 – 8/2020

**Project Title:** Plasmonic-dielectric metamaterials for enhanced spectroscopy and controlled chemistry.

**PIs:** N. Large (PI, 100%)

**Sponsor:** Office of Naval Research (ONR), U.S. Department of the Navy (DoN)

**Funded Amount:** $36,531

**Project Period:** 9/2017 – 9/2019
# Electron Energy-Loss Spectroscopy in Finite-Difference Time-Domain

## A. Related Publications

1. Cao, Y.; Manjavacas, A.; Large, N.; Nordlander, P. Electron energy-loss spectroscopy calculation in finite-difference time-domain package. *ACS Photonics* **2015**, *2*, 369–375.

2. Zhang, M.; Large, N.; Koh, A. L.; Cao, Y.; Manjavacas, A.; Sinclair, R.; Nordlander, P.; Wang, S. X. High-density 2D homo- and hetero-plasmonic dimers with universal sub-10-nm gaps. *ACS Nano* **2015**, *9*, 9331–9339.

3. Yazdi, S.; Daniel, J. R.; Large, N.; Schatz, G. C.; Boudreau, D.; Ringe, E. Reversible shape and plasmon tuning in hollow AgAu nanorods. *Nano Lett.* **2016**, *16*, 6939–6945.

4. Kraja, R.; Large, N. Computational modeling of the electron energy-loss spectroscopy (EELS) of periodic silver nanoparticle arrays. *Submitted* **2023**.

## B. Background

Energy-loss spectroscopy using fast electrons was used in the first experimental detection of surface plasmons in metals.⁵², ⁵³ Since these pioneering studies, electron energy-loss spectroscopy (EELS) has become a unique tool for probing surface plasmons of metallic nanostructures with unprecedented spatial ($<$ 1 nm) and energy ($<$ 100 meV) resolution. With the ability of electrons to probe dark surface plasmons and to provide quantitative information on surface plasmon kinetics, damping, and dispersion, EELS has become an irreplaceable tool for experimental study of surface plasmons. Consequently, EELS is widely used to investigate optical properties of complex metallic nanostructures (i.e., complicated geometries, strongly coupled nanosystems). The particular nature of the electromagnetic field of an electron (i.e., fast moving point charge, **Figure 1a**) makes the theoretical modeling of EELS data more involved than the modeling of optical experiments. For such simulations a myriad of different numerical techniques have been developed in the past years, including boundary element method (BEM),⁵⁴ discrete dipole approximation (DDA),⁵⁵ and finite element method (FEM). Although these methods are able to predict and interpret experimental EELS spectra, they possess some of the following significant drawbacks: (i) need for large computational resources, (ii) limitations to non-penetrating electron trajectories, (iii) requirement of highly symmetrical geometries, and (iv) complexity (i.e., programming skills required, absence of user friendly interface). These drawbacks are the main reason why EELS experimental results are generally compared to numerical simulations for optical excitations. The comprehensive review article by Kociak and Stéphan explicitly
states that, despite an implementation for cathodoluminescence (CL) spectroscopy in FDTD and the EELS implementation in the discontinuous Galerkin time-domain (DGTD) method, there is a clear lack of time-domain numerical methods for electron-beam spectroscopy calculations.⁵²

Here, I proposed a novel numerical procedure for EELS simulations employing a reliable and widely used commercial package: Ansys-Lumerical FDTD.⁵⁶ This package implements a high performance 3D solver that offers a user friendly environment for solving Maxwell's equations using the FDTD method. We benchmarked our method by direct comparison with results from the well-established BEM method⁵⁴ for three representative systems: (i) isolated nanospheres, (ii) nanosphere dimers, and (iii) nanodisks supported by a substrate. In the particular case of non-penetrating trajectories for nanospheres, we also compared our method with analytical results from Mie theory. Finally, in order to demonstrate the power of this implementation we performed EELS calculations for a more complex geometry consisting in a supported bowtie antenna for which experimental results are available.

### C. Method

The theoretical formalism for EELS simulations has been extensively discussed in the literature.⁵³ The physics of EELS can be understood as follows: the electric field produced by an electron moving with constant velocity $\mathbf{v}$ polarizes the nanostructure placed in the vicinity of its trajectory (cf. Figure 1a). This, in turn, induces an electric field $\mathbf{E}^{\text{ind}}$ that acts back on the electron, exerting a force that produces the energy loss, which can be written as

$$ \Delta E_{\text{EELS}} = e \int \mathbf{v} \cdot \mathbf{E}^{\text{ind}} [\mathbf{r}_e(t), t] dt = \int_0^\infty \hbar \omega \Gamma_{\text{EELS}}(\omega) d\omega, \quad (1) $$

where $e$ is the elementary charge, $\mathbf{r}_e(t)$ represents the electron trajectory, and

$$ \Gamma_{\text{EELS}}(\omega) = \frac{e}{\pi \hbar \omega} \int \text{Re}\{e^{-i\omega t} \mathbf{v} \cdot \mathbf{E}^{\text{ind}}[\mathbf{r}_e(t), t]\} dt + \Gamma_{\text{bulk}}(\omega) \quad (2) $$

is the energy-loss probability per unit of frequency $\omega$. The second term of this expression, $\Gamma_{\text{bulk}}$, represents the bulk loss probability. This contribution can be calculated using the following analytical expression

$$ \Gamma_{\text{bulk}}(\omega) = \frac{e^2 L}{\pi \hbar v^2} \text{Im} \left\{ \left( \frac{v^2}{c^2} - \frac{1}{\varepsilon} \right) \ln \left( \frac{c^2 q_c^2 / \omega^2 - \varepsilon}{c^2 / v^2 - \varepsilon} \right) \right\}, \quad (3) $$

where $L$ is the length of the electron trajectory inside the medium, $\varepsilon(\mathbf{r}, \omega)$ is the medium dielectric permittivity, and $q_c \approx \hbar^{-1}[(m_e v \phi_{\text{out}})^2 + (\hbar\omega/v)^2]^{1/2}$ is the cutoff momentum determined by the electron mass $m_e$, the electron velocity $v$, and the collection angle $\phi_{\text{out}}$ of the microscope. This expression is valid within the local response approximation, in which only low enough momentum transfers below $q_c$ are collected. In the remainder of this chapter, we do not consider the bulk contribution. In addition, and without loss of generality, we assume the electron trajectory to be in the $(x, z)$ plane, parallel to z-axis, and separated from the origin by the impact parameter $b$, so $\mathbf{v} = v\hat{\mathbf{z}}$ and $\mathbf{r}_e(t) = (b, 0, vt)$.

The electric field created by an arbitrary current distribution (such as a beam of electrons) can be written as

$$ \mathbf{E}[\mathbf{r}(t), \omega] = -4\pi i\omega \int \mathbf{G}(\mathbf{r}, \mathbf{r}', \omega)\mathbf{j}(\mathbf{r}', \omega)d\mathbf{r}', \quad (4) $$

in terms of the Green tensor of Maxwell's equations $\mathbf{G}(\mathbf{r}, \mathbf{r}', \omega)$ (note that we use Gaussian units). This quantity is defined as the solution of Helmholtz equation:
$$ \left[-\nabla \times \nabla \times + \varepsilon(\mathbf{r}, \omega)\frac{\omega^2}{c^2}\right] \mathbf{G}(\mathbf{r}, \mathbf{r}', \omega) = \frac{1}{c^2}\delta(\mathbf{r} - \mathbf{r}')\mathbf{I}, \qquad (5) $$

where **I** is the unit tensor. Using the electron current density relevant for EELS, $j(\mathbf{r}', \omega) = -e\delta(x-b)\delta(y)e^{i\omega z/v}\hat{\mathbf{z}}$, and Eq. (3) we can rewrite the loss probability [Eq. (2)] as

$$ \Gamma_{\text{EELS}}(\omega) = \frac{4e^2}{\hbar} \iint_{-\infty}^{\infty} \cos\left[\frac{\omega(z-z')}{v}\right] \operatorname{Im}[-G_{zz}^{\text{ind}}(z,z',\omega)]dzdz'. \qquad (6) $$

Here, $G_{zz}^{\text{ind}}(z, z', \omega) = \hat{\mathbf{z}} \cdot [\mathbf{G}(\mathbf{r}, \mathbf{r}', \omega) - \mathbf{G}^0(\mathbf{r}, \mathbf{r}', \omega)] \cdot \hat{\mathbf{z}}$, with $\mathbf{G}^0(\mathbf{r}, \mathbf{r}', \omega)$ being the Green tensor for vacuum (to simplify the notation we have omitted the lateral spatial coordinates $x=b$ and $y=0$ in the arguments of the Green functions and related quantities). In the derivation of this expression, we have also used the reciprocity property of the Green tensor $\mathbf{G}(\mathbf{r}, \mathbf{r}', \omega) = \mathbf{G}^\top(\mathbf{r}', \mathbf{r}, \omega)$. Interestingly, $G_{zz}^{\text{ind}}(z, z', \omega)$ can be obtained from the z-component of the electric field induced at position z, by an electric dipole of amplitude $p(z', \omega)$ placed at $z'$ and oriented along z-axis

$$ G_{zz}^{\text{ind}}(z, z', \omega) = -\frac{1}{4\pi\omega} \frac{E_z^{\text{ind}}(z, \omega)}{p(z', \omega)} \qquad (7) $$

Using this expression, we can rewrite Eq. (6) as

$$ \Gamma_{\text{EELS}}(\omega) = \frac{e^2}{\pi \hbar \omega} \iint_{-\infty}^{\infty} \cos\left[\frac{\omega(z-z')}{v}\right] \operatorname{Im}\left[\frac{E_z^{\text{ind}}(z, \omega)}{p(z', \omega)}\right] dz dz' \qquad (8) $$

Therefore, all we need to do in order to obtain the loss probability is to compute the induced electric field generated by an electric dipole along the electron trajectory. Notice that Eq. (4) involves only the imaginary part of the induced field, which remains finite even at the position of the dipole. The use of the induced field, in place of the total field, for non-penetrating trajectories is not required by the theoretical formalism, since an electron cannot produce energy loss in absence of material structures. However, due to the finite accuracy of the FDTD calculations we choose to work with the induced field in order to minimize any numerical instability originating from the calculation of the fields. In order to calculate the field, we can employ any Maxwell's equation solver. Here, we choose to work with the commercial FDTD package from Ansys-Lumerical due to its convenient user environment.⁵⁶ The computation procedure starts by setting a 3D FDTD simulation domain with perfectly matched layers (PMLs) to prevent spurious reflections from outer boundaries (Figure 1b,d). We then insert the nanostructure and define an override mesh that allows us to manually adjust the mesh grid size in a particular region. This allows us to optimize the discretization of the physical object and to improve the convergence. After that, we place a 1D (linear) monitor along the electron trajectory that allows us to calculate the electric field at specific points. The monitor is extended across the entire simulation domain through the PMLs. Again, to improve the convergence and to ensure a proper spatial discretization along the electron trajectory, we place a second override mesh on top of the monitor. Next, we position an electric point dipole $p(z', \omega)$ on the mesh grid points, aligned with the electron path, and successively displace it from mesh point to mesh point along the electron trajectory from $z' = z_{\text{max}}$ (upper PML) to $z' = z_{\text{min}}$ (lower PML). This electric dipole acts as a source in Maxwell's equations. For each position, we record the z-component of the total electric field $E_z(z, \omega)$ along the entire monitor (Figure 1b). To obtain the induced electric field $E^{\text{ind}}(z, \omega)$ we subtract the background electric field $E_z^0(z, \omega)$ from $E_z(z, \omega)$. The former quantity is calculated using the previously described protocol but removing all the physical objects (e.g., nano-object and substrate) from the simulation domain (Figure 1c). For situations where the electron trajectory penetrates into the absorbing medium (e.g., the metal), one needs to
calculate the corresponding background electric field generated by an electric dipole placed
in an infinite space filled with the corresponding material. In practice, this can be
accomplished by placing the dipole at the center of a sphere composed of this material,
whose diameter must be chosen large enough to minimize the field spill out into the
surrounding medium. Here we choose this diameter equal to the monitor length. In metals, as
the field of the electric dipole decays to zero after a few tens of nanometers, no field exits the
micron-sized metallic domain, which thus can be considered as infinite from the dipole's
point of view. Incidentally, since FDTD simulations are not stable when an absorbing
medium is extended through the PMLs, this forces us to enlarge the simulation domain and
the monitor to prevent the absorbing medium to reach into the PMLs. Once the induced
electric field is calculated, the EELS spectrum is readily computed using Eq. (8). Finally we
note that the velocity of the electron only enters in Eq. (8) through the cosine function. This
allows us to compute EELS spectra for any electron velocity from a single FDTD calculation.
It is important to notice that one can either calculate all the dipole positions in one single
FDTD calculation or split each dipole position into smaller sub-calculations to increase the
parallelization and optimization. The results presented in this chapter are performed using the
later procedure.

Lastly, it is worth mentioning that, unless DDA and BEM, the proposed FDTD method does
not use a direct implementation of the electron source. While DDA and BEM use a rigorous
description of the electric field generated by a swift electron as the source term in the
governing equations,⁵⁵,⁵⁷ EELS-FDTD uses a Green’s tensor approach and makes use of
electric point dipoles as sources. In other words, current EELS implementations in BEM and
DDA do not use Eq. (8).

D. Results and Discussion

1. *Isolated Nanostructure*

We first illustrate the EELS-FDTD implementation for the case of an isolated gold
nanosphere of diameter *a* = 160 nm placed in vacuum. To benchmark our method, BEM
calculations are performed using an electron source implemented in the axial symmetry
version of this semi-analytical method following the formalism established by García de
Abajo and Aizpurua.⁵⁷ In addition, the simple spherical geometry allows us to perform
analytical (Mie theory). We use the dielectric function of gold tabulated by Johnson and
Christy.⁵⁸ While the experimental values are used as is in the BEM and Mie calculations,
analytical multi-coefficient models (MCMs) are used in FDTD to fit these experimental data
and overcome the difficulty of adapting spectrally tabulated dielectric permittivities into
time-domain methods. The nanoparticle center is placed at the origin of the coordinate
system. To ensure a good convergence, we set a monitor length of 1500 nm (i.e., z<sub>max(min)</sub> =
±750 nm). The other parameters used in the Ansys-Lumerical FDTD simulation are set as
follows: a simulation time of 100 fs with an auto shutoff parameter of 10<sup>-5</sup>, a mesh accuracy
of 5 (i.e., 22 mesh points per wavelength), and mesh refinement algorithm set to "conformal
variant 0" allowing for a non-uniform mesh over the FDTD domain. An initial simulation is
performed to calculate the total electric field E<sub>z</sub> at each point along the electron trajectory in
the frequency range 1–3 eV. A second calculation is then performed for the same electric
dipole positions in absence of the nanoparticle to calculate the background electric field
generated by the dipole in vacuum E<sub>z</sub><sup>0</sup>. Schematics of these two configurations are shown in
**Figure 1b,c**. Then, the induced electric field is computed as E<sub>z</sub><sup>ind</sup>(z, ω) = E<sub>z</sub>(z, ω) –
E<sub>z</sub><sup>0</sup>(z, ω) and inserted into Eq. (8) to obtain the EELS spectra. The results of this calculation
are shown with blue lines in **Figure 2** for three different impact parameters: *b* = 120 nm
(away from the nanoparticle, bottom), *b* = 82 nm (in close proximity to the nanoparticle,
center), and $b = 0$ nm (through the center of the nanoparticle, top). In all the cases we use an electron speed equal to half of the speed of light (i.e., $v = 0.5c$), which corresponds to a kinetic energy of 80 keV.

<figure><img src="image_13.png" /><figcaption>**Figure 1.** (a) Schematics showing the electric field of an electron moving with velocity $v$ along a straight line trajectory separated from the nanostructure by an impact parameter $b$. (b,c) Geometry used to calculate the induced electric field $E^{\text{ind}} = E_z - E_z^0$. Here, $E_z$ and $E_z^0$ are the $z$-component of the electric field generated by an electric dipole $p(z', \omega)$ at position $z$ in presence and in absence of the metallic nano-object, respectively. (d) FDTD simulation setup showing: the simulation domain with PML boundary conditions, the nanostructure, the electric dipole source, the 1D monitor used to record the electric field along the electron path, and the override meshes used to discretize the nanostructure and the electron trajectory.</figcaption></figure>

The results obtained with the EELS-FDTD implementation (blue lines) show a strong peak at 2.4 eV, in very good agreement with BEM (red lines) and Mie theory (black triangles) calculations. This peak corresponds to the quadrupolar mode of the nanoparticle. The dipolar mode of the nanoparticle only appears as a shoulder in the spectrum at smaller energies. Interestingly, the position of this peak depends on the impact parameter due to retardation effects, originating from the frequency dependence of the field from the electron. In contrast to several published methods,⁵⁴,⁵⁵ the EELS-FDTD implementation can also handle penetrating trajectories (Figure 2, top). However, for such cases one has to be careful when performing the FDTD simulation.
<figure><img src="image_14.png" /><figcaption>**Figure 2.** EELS spectra for a gold nanosphere of diameter *a* = 160 nm calculated with the EELS-FDTD implementation (blue lines) for three impact parameters. From top to bottom: *b* = 0 nm, *b* = 82 nm, and *b* = 120 nm. The EELS-FDTD calculations are compared with BEM calculations (red lines), and with Mie theory (black triangles, only for *b* > *a*/2). The spectra for *b* = 120 nm are multiplied by 10 to improve the clarity of the figure. The electron speed is 0.5c (i.e., 80 keV) in all cases.</figcaption></figure>

The relative position of the electric dipoles with respect to the nanoparticle surface can artificially introduce numerical errors. We show that when an electric dipole is placed exactly at the nanoparticle surface, it produces an overestimation of the EELS signal. This is easily solved by slightly displacing the entire nanostructure along the *z*-axis with respect to the monitor mesh grid (typically 1/3 mesh step). The discrepancies observed between FDTD and BEM results (Figure 2, top) can be minimized by increasing the number of dipoles used in the calculation, in particular in the part of the trajectory close to the nanostructure.

Although our EELS-FDTD implementation requires performing of a large number of short sub-calculations (*i.e.*, one per dipole position), it allows for reaching a better convergence level at lower computational cost than DDA in specific configurations. BEM, due to the axial symmetry nature of this particular problem, allows the user to perform the same calculation much faster. A comparison of the computational resource used by FDTD, DDA, and BEM to calculate the EELS spectrum for an impact parameter *b* = 82 nm (Figure 2, center) is provided in Table 1.
<table><caption>**Table 1.** Comparison of the computational resource (time, memory) between FDTD, DDA, and BEM. The EELS spectra are calculated for an isolated gold nanosphere of diameter *a* = 160 nm with an impact parameter *b* = 82 nm, for an energy range from 1 to 4 eV, and an energy resolution of 15 meV. The discretization parameter *d* corresponds to the minimum volume mesh size (FDTD), the boundary element size (BEM), and the interdipole distance (DDA). The total time is calculated for each EELS calculation when the sequences are run in serial.</caption><thead><tr><th></th><th><i>d</i><br>[nm]</th><th>Number of<br>elements</th><th>Number of<br>sequences</th><th>Average time<br>per sequence<br>[min]</th><th>Total<br>compute time<br>[h]</th><th>Average<br>memory<br>[MB]</th></tr></thead><tbody><tr><td>FDTD #1</td><td>1</td><td>6,011,499<br>Yee's cells</td><td>750 dipole<br>positions</td><td>8.6 ± 2.1</td><td>~ 105</td><td>2,874</td></tr><tr><td>FDTD #2</td><td>1</td><td>1,642,560<br>Yee's cells</td><td>500 dipole<br>positions</td><td>4.6 ± 0.6</td><td>~ 35</td><td>1,659</td></tr><tr><td>DDA</td><td>0.8</td><td>4,188,896<br>dipoles</td><td>200 energy<br>values</td><td>62 ± 20</td><td>> 210</td><td>15,135</td></tr><tr><td>BEM</td><td>2.5</td><td>100 boundary<br>elements</td><td>200 energy<br>values</td><td>6.8 s</td><td>~ 25 min</td><td>14.50</td></tr></tbody></table>

## 2. Interacting Nano-Objects

Nanoparticle dimers have been extensively studied using EELS.⁵² For this reason they constitute another ideal system to benchmark our EELS-FDTD method. Here, we study a dimer of closely spaced (i.e., strongly interacting) silver nanospheres of diameter *a* = 160 nm placed in vacuum. The gap size is fixed to *g* = 5 nm, and we use the dielectric function for silver tabulated by Palik⁵⁹ in the BEM calculations and MCMs fit of the later in the FDTD calculations. The other simulation parameters are the same as for the gold nanosphere.

**Figure 3** shows the results obtained with our EELS-FDTD method (blue lines) for three different impact parameters: *b* = 164.5 nm (dimer end, top), *b* = 82.5 nm (through the center of one of the NP, center), and *b* = 0 nm (in the gap, bottom). In all the cases we use an electron speed of *v* = 0.5c (i.e., 80 keV). It is well known that for large gaps the nanoparticles interact only weakly and the resulting dimer plasmons are essentially bonding and anti-bonding combinations of the nanoparticle plasmons of the same multipole order (e.g., *l* = 1, dipole). Here, due to the very small gap-to-diameter ratio (*g*/*a* = 0.03) the plasmon modes of the dimer contain contributions from various multipole orders. When the electrons pass in close proximity to the dimer end (upper spectra), both FDTD and BEM show the appearance of two weak localized surface plasmon resonances (LSPRs) at 1.5 and 2.3 eV and stronger features at 3.25 and 3.5 eV (i, ii, iii, and iv, respectively). Though mixed with higher values of *l* the 3.5 eV is dominated by the anti-bonding dipole while the 3.25 eV feature results from the strong hybridization of high order modes. The low energy and weak features observed at 1.5 and 2.3 eV are the longitudinal bonding dipole and bonding quadrupolar modes, respectively. The EELS maps associated to the four LSPRs (i-iv) are shown in **Figure 3b** with a 4 nm spatial resolution. Here, for the sake of simplicity and for computational considerations, we excluded penetrating trajectories. Although the present dimer is much larger (i.e., more retardation effects) and has a smaller gap-to-diameter ratio (i.e., more hybridization), the results are in good agreement with EELS measurements reported in literature. The large energy splitting between the bonding and anti-bonding dipolar LSPRs is the signature of a strong coupling. When the electron beam follows a penetrating trajectory (center spectra) the high energy anti-bonding dipole (3.5 eV) remains strong and the bonding
quadrupole (2.3 eV) is strengthened. For a penetrating electron trajectory, also the bulk plasmon mode at 3.8 eV can be excited. However, due to the present numerical procedure to calculate the induced electric field the bulk plasmon cannot be observed in our EELS-FDTD method. This procedure does not affect the localized surface plasmon modes, which are determined by the surface charges around the structure. For the gap center trajectory (lower spectra) the bonding dipolar dimer mode cannot be excited and the spectrum is dominated by the anti-bonding dipolar LSPR. Incidentally, due to retardation effects the position of this mode is slightly different for the different excitation configurations.

<figure><img src="image_15.png" /><figcaption>**Figure 3.** (a) EELS spectra for a dimer of silver nanospheres of diameter *a* = 160 nm separated by a gap *g* = 5 nm. Results obtained with the EELS-FDTD implementation (blue lines) are compared with BEM calculations (red lines) for three different impact parameters: from top to bottom: *b* = 164.5 nm, *b* = 82.5 nm, and *b* = 0 nm. In all cases the electron speed is fixed to 0.5c (i.e., 80 keV). Inset: Zoom-in view of the low-energy part of the upper spectrum. (b) EELS maps calculated at four different energies.</figcaption></figure>

### 3. Supported Nanostructures

EELS experiments require very thin non-absorbing substrates to minimize energy losses. Typical EELS substrates are made of mica, silica ($SiO_2$), silicon nitride ($Si_3N_4$, $SiN_x$), or carbon (C), and their thickness generally ranges from 5 to 50 nm. Even though such thin substrates introduce a negligible EELS background and are often considered to have a negligible effect on the optical properties of the supported nanostructure (i.e., small spectral shift with respect to the free-standing nanostructure), they can be crucial in some situations.
For this reason, we calculate the EELS spectra for gold nanodisks of diameter *a* = 50 nm and height *h* = 15 nm placed on 30 nm thick substrates of different dielectric permittivities: ε = 1 (free standing); ε = 2 (SiO₂); and ε = 4 (Si₃N₄). As in the case of the isolated sphere, we use the gold dielectric function tabulated by Johnson and Christy⁵⁸ in the BEM calculations and MCMs fit of the later in the FDTD calculations. The FDTD simulation parameters are the same as in previous cases and the electron speed is set to v = 0.5c (80 keV). The results are shown in **Figure 4** for FDTD (blue lines) and BEM (red lines) for an impact parameter *b* = 27 nm (in close proximity to the nanodisk). As expected, the dipolar LSPR (2.3 eV for ε = 1) red-shifts with increasing permittivity of the substrate (1.97 eV for ε = 4). Interestingly, it has to be noticed a change in the LSPR lineshape for ε = 4. This effect, along with the spectral shift, clearly shows that even very thin substrates can have a significant impact on the EELS spectrum and may need to be included in the simulations.

Finally, to highlight the power and the flexibility of our EELS-FDTD implementation, we perform EELS calculations for a supported bowtie antenna. This complex structure is composed of two gold equilateral triangles with a lateral length *a* = 80 nm, a height *h* = 15 nm, separated by a gap *g* = 4 nm. The gold bowtie structure is placed on top of a 30-nm thick SiN substrate. We also included a 2.5-nm thick chromium (Cr) adhesion layer. The dielectric permittivity of SiN is assumed to be constant and equal to 5.5, while the corresponding one for Cr is described by MCMs fit of tabulated data.⁵⁸

<figure><img src="image_16.png" /><figcaption>**Figure 4.** EELS spectra for a gold nanodisk of diameter *a* = 50 nm and height *h* = 15 nm placed on top of a 30 nm thick substrate with dielectric permittivity ε = 1, 2, and 4 calculated with the EELS-FDTD implementation (blue lines), and BEM (red lines). The impact parameter is fixed to *b* = 27 nm and the electron speed is 0.5c (i.e., 80 keV).</figcaption></figure>

The results of this simulation are shown in **Figure 5**. There, we observe that for edge excitations (**Figure 5a**) the spectra display two LSPRs (i, iv) located at 1.27 and 2.39 eV, respectively. In this case, different impact parameters produce very different intensity ratios for these LSPRs. On the other hand when electron trajectory crosses the center of the gap (**Figure 5b**) the EELS spectrum displays two distinct LSPRs (ii, iii) at 1.68 and 2.17 eV, respectively, with relative intensities that are much less dependent on the impact parameter. When displacing the electron trajectory off-axis, the LSPR at 2.17 eV (iii) becomes weaker progressively, while the one at 1.68 eV (ii) remains unchanged. All these results, along with
the results for a single triangular prism, are in good agreement with the experimental observations by Yang and co-workers.⁶⁰,⁶¹ The small discrepancies, mainly in the linewidths, are related to the presence of a thicker Cr adhesion layer in our calculations, which is known to introduce a broadening and a red-shift of the LSPRs.⁴⁶,⁴⁷

<figure><img src="image_17.png" /><figcaption>**Figure 5.** EELS spectra for a gold bowtie antenna calculated with the EELS-FDTD implementation. Each triangle has a lateral length *a* = 80 nm, height *h* = 15 nm, and gap *g* = 4 nm. The bowtie antenna is supported by a 30 nm thick SiN substrate. A 2.5 nm Cr adhesion layer is included. The electron speed is taken equal to 0.5c (i.e., 80 keV) in all cases. (a) EELS spectra for edge excitation with four different impact parameters: 2 (blue), 5 (purple), 10 (red), and 15 nm (black) (cf. inset). (b) EELS spectra for gap excitation with five different impact parameters: 0 (blue), 2 (purple), 5 (red), 10 (brown), and 15 nm (black) (cf. inset). (c) EELS maps calculated at (i) 1.27, (ii) 1.68, (iii) 2.17, and (iv) 2.39 eV.</figcaption></figure>

**Figure 5c** shows the EELS maps corresponding to the four LSPRs (i-iv) calculated using our EELS-FDTD implementation with a 2 nm spatial resolution. Similarly to the dimer, we here choose to exclude penetrating trajectories. The nature of the LSP modes can straightforwardly be determined from these maps. Modes (i) and (ii) correspond to the dipolar bonding (bright) and anti-bonding (dark) modes, respectively. The map for 2.39 eV shows strong loss probability from each bowtie edge underlining the high order nature of mode (iv). These results are also in excellent quantitative agreement with the recent studies by Yang and co-workers.⁶⁰,⁶¹ Finally, mode (iii) displays a loss probability which is strongly localized at the bowtie gap. This spatial confinement directly correlates with the rapid vanishing of mode (iii) when the electron trajectory is displaced off gap (Figure 5b).

## 4. Complex Nanostructures

The new EELS-FDTD method was later successfully applied to more complex, realistic nanostructures. The first complex system of interest consisted in high-density homo- (Au-Au and Ag-Ag) and hetero- (Au-Ag) dimer substrates with sub-10-nm gaps. Dr. Shan Wang from Stanford University fabricated these high-density, large-scale nanostructured substrates using a novel approach based on nanoimprint lithography and shadow evaporation.⁶² The
dimers were modeled as two slightly different hemi-ellipsoids of center height $h = 40$ nm, with a gap $g = 5$ nm, on top of a $Si_3N_4$ membrane. The in-plane dimensions of the two ellipsoids are $L_1, L_2$ for the long axes and $l_1, l_2$ for the short axes, respectively.

For the Au-Au dimer, the long/short axis lengths for the two nanoparticles are 190/150 nm and 170/130 nm, respectively (Figure 6a,b). For the Ag-Ag dimer, these values were set at 160/130 nm and 155/110 nm, respectively (Figure 6c,d). The calculated EELS spectra shown in Figure 6 (dashed lines) are in good agreement with the experimental results and have captured the salient features observed in the EELS measurements (solid lines). More specifically, EELS-FDTD calculations on the Au-Au dimer reveal LSPRs around 1.2 and 1.6 eV for the edge excitation (Figure 6a) and several LSPRs between 1.6 and 2.1 eV for the gap excitation (Figure 6b). Similarly, the calculations for the Ag-Ag dimer reveal LSPRs around 1.3, 1.8–2.1, and 3.3–3.6 eV for the edge excitation (Figure 6c) and from 1.8 to 2.2 eV, and at 3.6 eV for the gap excitation (Figure 6d). It is important to notice that the EELS-FDTD spectra of the Au-Au and Ag-Ag dimers exhibit multiple smaller peaks with spectral separations approaching the experimental resolution (150 meV), which result in the broad bands observed experimentally.

As the dimers are formed from elongated nanoparticles, the position of the LSPRs depends on the geometrical aspect ratios. Variations in aspect ratios across the entire sample also contribute to inhomogeneous broadening of the LSPR when collectively probed by optical spectroscopy. However, EELS is much less sensitive to this source of inhomogeneous broadening as it probes individual nanoparticles. Other morphological details such as surface roughness, small height inhomogeneity, and deviation from the ideal semi-ellipsoidal geometry used in the theoretical EELS-FDTD calculations contribute to LSPR broadening and LSPR shift. Contrary to ensemble measurements, fine morphological details play a significant role in EELS of individual nanostructures due to the highly localized nature of the electron beam. Local irregularities (e.g., curvature and indentations) in the gap region strongly influence the local density of states in this particular region and therefore, are reflected in the associated EELS spectrum. Also, high-angle annular dark-field (HAADF) images revealed that tiny nanometer-size dots were formed around the dimers. The change in surrounding environment also introduces spectral shifts of the observed LSPRs. Such morphological features can be accounted for in the EELS-FDTD calculations by introducing very fine sub-nanometer FDTD mesh. However, we decided to keep the geometry ideally smooth, thus allowing us to investigate the physics of the system while keeping the computational cost reasonable.

Similar analysis was also conducted on the Au-Ag hetero-dimer with size parameters ($L_{Au}/l_{Au}$, $L_{Ag}/l_{Ag}$) = (175/100, 175/110) nm (Figure 6e). The experimental and calculated EELS spectra are shown in Figure 6f. The EELS spectrum from the Au end (red) exhibits two plasmon bands between 1.5 and 2 eV while the spectrum from the Ag end (green) exhibits two plasmon bands between 1.5 and 2 eV and smaller LSPRs around 3.5 eV. The EELS spectrum from the gap center (blue) is dominated by the Ag component of the heterodimer, exhibiting a broad LSPR around 2 eV and a weaker LSPR at 3.5 eV.

Another complex system of interest where EELS-FDTD was successful at interpreting experimental EELS results consists of hollow AgAu-alloy nanorods created by partial galvanic replacement. The fabrication and experimental EELS characterization of these nanostructures was done by Dr. Emilie Ringe from the University of Cambridge.⁶³ The modeled geometry and dimensions were based on STEM-HAADF images of the nanorods (Figure 7a), and were modeled as cylinders with semi-spherical ends, with a length $L_{Au} = 195$ nm and a radial diameter $a_{Au} = 47$ nm. The radially concentric water void was also
modeled cylindrical with semi-spherical ends with $L_{water} = 80$ nm and $a_{water} = 35$ nm. Its relative longitudinal position with respect to the center of the gold nanorod is given by $\Delta$. While the complete study looked at various $\Delta$ values,⁶³ **Figure 7** only shows the results for $\Delta = 0$ nm (centered void).

<figure><img src="image_18.png" /></figure>

**Figure 6.** (a,b) EELS of Au-Au dimer with electron beam directed (a) on the edge and (b) in the gap of the dimer. The solid lines represent the experimental data and the dashed lines represent the theoretical results. (c,d) EELS of Ag-Ag dimer with electron beam directed (c) on the edge and (d) in the gap of the dimer. The solid lines represent the experimental data and the dashed lines represent the theoretical results. The insets are the corresponding HAADF images with the red dots indicating the position of the electron beam. (e) HAADF image of a Au-Ag heterodimer, with the left part made of Au and the right part made of Ag. (f) EELS spectra (solid lines) and EELS-FDTD spectra (dashed lines) were acquired with electron beam incident on the Au edge (red), Ag edge (green), and gap (blue).

The calculated EELS spectra (not shown here) were in excellent agreement with the experimental EELS measurements. Decomposition of the calculated EELS spectrum revealed the presence of various plasmon modes between 1.27 and 2.37 eV. The EELS-FDTD loss probability maps (**Figure 7b**) together with surface charge density calculations (**Figure 7c**) allowed for revealing the nature of the different plasmon modes. The calculated loss probability distributions were in excellent agreement with the experimental EELS results (**Figure 7a**).
<figure><img src="image_19.png" /><figcaption>**Figure 7.** (a) STEM-HAADF images of a hollow AuAg nanorod and experimental EELS maps of each of the four LSPRs. (b) Decomposition of the calculated loss probability [$\Gamma_{EELS}$, Eq. (8)] maps calculated using the EELS-FDTD method for the center state nanorod (i.e., $\Delta = 0$ nm). (c) Calculated surface charge density distribution (red: positive charge, blue: negative charge) for each of the four LSP modes obtained in the FDTD calculations from panel b.</figcaption></figure>

Specifically, the modes at 2.1 and 2.37 eV in EELS-FDTD appear to be the bonding and antibonding combinations of the transverse dipoles, respectively. These two modes, although spectrally shifted from the experimental values, correspond to the modes experimentally observed at 1.8 and 2.5 eV, respectively. Additionally, a quadrupolar mode can be seen at 2.28 eV (2.1 eV experimentally). Finally, the decomposition of the EELS-FDTD spectrum contains two modes close in energy (1.27 and 1.33 eV). These energies are well aligned with the 1.2 eV mode observed experimentally and could be attributed a splitting of the bonding and antibonding combinations of the longitudinal dipoles, respectively. Further, mapping the nanorod LSPRs in EELS has revealed a strong interaction between the void and nanorod localized surface plasmons.

## 5. Periodic Nanoparticle Arrays

A recent theoretical work by Zundel and Manjavacas showed that under an optical excitation, a 5×5 NP array exhibit the same optical response as a purely infinite array.⁶⁴ While this result may appear consistent with the collective excitation of the finite array, localized excitations may induce a different response of the array. Here, we performed EELS calculations to investigate the finite-size effects of periodic arrays of Ag nanoparticles. Although EELS has also been successfully used to study photonic crystals³² and has become more popular for single NPs or small NP assemblies,⁵³ EELS studies of 2D and 3D periodic systems remain very scarce.³², ⁶⁵, ⁶⁶ We calculated the loss probability spectra and associated EELS maps of 2D periodic arrays with dimensions ranging from a small 2×2 to a 10×10 nanoparticles. The nanoparticles were modeled as spheres of 40 nm in diameter and with 2 nm gaps (Figure 8).

Because EELS allows for the excitation of optically dark modes, the EELS spectra exhibit larger number of LSP modes as compared to the optical extinction spectrum (Figure 8, left). When computing the loss probability and surface charge distributions for selected LSPRs (Figure 8, right), we can clearly see that these resonances can be associated with edge and corner modes. In Figure 8i-iii, for instance, both the EELS maps and surface charge distributions clearly show a strong localization of the LSP modes along the edges of the array. On the other hand, the mode at 2.55 eV (Figure 8iv) is more strongly localized at the
corners of the array. The strong localization of these LSP modes along the edges and corners can only occur in a finite-size array. Similar behavior was observed for other array sizes. These finite-size effects can be observed only because EELS make use of a local excitation, which allows for selectively exciting individual modes, including optically dark modes. Our work revealed that contrary to what Zundel and Manjavacas concluded in their recent work, finite-size arrays larger than 5×5 do not behave as infinite arrays, but rather sustain a larger variety of edge and corner modes which can be probed in electron-based spectroscopies (EELS, CL) and can play an important role in applications such as sensing.

<figure><img src="image_20.png" /><figcaption>**Figure 8.** **Left:** Loss probability spectra calculated for a 6×6 NP array. The EELS spectra are color-coded with the various impact parameters shown in the inset. The extinction cross-section calculated for a plane wave excitation is shown as a comparison (blue). **Right:** Loss probability maps and surface charge density distributions calculated at 4 different energies.</figcaption></figure>

## E. Conclusions

We have presented a simple procedure to calculate the energy-loss probability of fast electrons interacting with metallic nanostructures. Although this method can be implemented with any Maxwell's equation solver we have chosen here to work with the commercial package Ansys-Lumerical FDTD due to the flexibility of the FDTD method, its user friendly environment, and its wide use. Contrary to most of the well-established methods, we have shown that this implementation can deal with both penetrating and non-penetrating trajectories and nanostructures of arbitrary geometries and morphologies, including substrates and adhesion layers. We have benchmarked our EELS-FDTD implementation by comparing the results with the well-established BEM method for different representative nanostructures such as nanospheres, nanoparticle dimers, and a nanodisk supported by a substrate. Furthermore, we have applied this method to study the EELS properties (spectra and loss probability distributions) of complex systems consisting of supported bowtie antennas, highly coupled homo- and hetero-dimers, and periodic nanoparticle arrays. Our EELS-FDTD method provides a simple and convenient approach for the calculation of EELS spectra and map from complex nanostructures of arbitrary shape and composition.
# Photonic Band Structures of Finite Periodic Nanostructures

## A. Related Publications

1. Montaño-Priede, J. L.; Large, N. Photonic band structure calculation of 3D finite nanostructured supercrystals. *Nanoscale Adv.* **2022**, *4*, 4589–4596.

2. Kraja, R.; Large, N. Computational modeling of the electron energy-loss spectroscopy (EELS) of periodic silver nanoparticle arrays. *Submitted 2023*.

## B. Background

Plasmonic materials exhibit limiting physical effects associated with their materials composition, such as inherent losses.⁶⁷⁻⁶⁹ A route to overcoming these limitations and achieving sustainable materials with improved optical properties consists of structuring the plasmonic building blocks in a well-organized manner,⁷⁰ forming one dimensional (1D) chains,²⁶⁻²⁸ two dimensional (2D) arrays,²⁸, ²⁹, ⁷¹⁻⁷³ or three dimensional (3D) superlattices.²⁸, ³⁰, ³¹, ³⁶, ⁷⁴⁻⁸¹ These supercrystals possess the unique ability to retain the intrinsic properties of their building blocks but also display unique collective properties originating from interparticle coupling effects. In recent years, DNA-linked 3D superlattices composed of metallic NPs have been investigated.³⁰, ⁷⁴, ⁷⁶, ⁷⁷, ⁸¹⁻⁸⁴ In their pioneering work, Mirkin and Schatz have shown that such 3D supercrystals exhibit strong novel optical behaviors resulting from the coexistence of plasmonic and photonic modes.²⁸, ³⁵, ⁷⁴, ⁸⁵ By continuously varying the NPs orientation, periodicity, spacing, size, and shape, as well as the crystal habit, an enormous set of fundamentally interesting metamaterials can be designed. With the countless superlattices now synthetically realizable,³¹ computational methods and theoretical models play a crucial role in identifying the supercrystals that exhibit the most exciting properties (Figure 9).³⁶, ⁸⁶ When the supercrystal lattice parameter is much smaller than the wavelength, the medium usually behaves like a uniform, continuous material, characterized by an effective dielectric permittivity and magnetic permeability, which are modelled using an effective medium theory such as the Maxwell-Garnett⁸⁷, ⁸⁸ and Bruggeman⁸⁹ models. However, when the lattice parameter is comparable to the optical wavelength, the photonic band structure plays a key role. In such assemblies, the interparticle distance becomes the critical parameter that governs the collective optical properties. When the distance is larger than the wavelength of optical excitation the NPs are uncoupled and the collective optical response is the incoherent summation of individual responses of each NP constituting the assemblies. Additional photonic behaviors arise when the interparticle distance satisfies Bragg's law (e.g., Rayleigh anomalies) or when multipolar interactions are involved. Finally, geometrical aspects such as the dimensionality and the habit of the supercrystal, the lattice parameter in each spatial direction, and the morphology of the NPs composing the
supercrystal play a key role in defining their optical properties including (an)isotropy/birefringence,⁸⁵,⁸⁶ the emergence of Fabry-Pérot resonances and whispering gallery modes (WGM),³⁶ and decrease of the photonic mode lifetime.

<figure><img src="image_21.png" /><figcaption>**Figure 9.** Approaches to calculate the properties of photonic crystals, i.e., infinite lattice for band structure calculations, effective medium theory for optical properties like reflectance, and proposed FDTD simulation for finite supercrystals.</figcaption></figure>

Understanding the optical properties of realistic, experimentally fabricated nanostructures require a deep understanding of morphological defects,⁹⁰⁻⁹² size distributions,⁹³ materials inhomogeneity, size effects,⁹⁴ and finite-size effects.⁶⁴, ⁹⁵, ⁹⁶ Computationally, periodic systems are very often regarded as infinite *via* the use of periodic boundary conditions applied to the simulation domain⁴², ⁸⁵, ⁹⁷ or effective medium approaches.³⁶, ⁷⁴, ⁹⁸ Only a few studies have focused on finite-size effects in 1D linear chains of nanoparticles,⁹⁹ and 2D periodic arrays.⁶⁴, ⁹⁶ All of these studies focused on very specific optical properties such as scattering, absorption, transmission, or reflection. With the recent advances in nanoscale and microscale fabrication techniques, plasmonic superlattices and their optical properties have rapidly become an interesting platform of novel optical properties. While semiconductor superlattices have been well studied and theoretical tools have been developed to study their optical properties,¹⁰⁰ the field of Plasmonics lacks computational models and methods to investigate these finite-size effects.¹⁰¹ On one hand, the nanoscale nature of the building blocks requires very fine spatial discretization of the computation domain to describe the near-field nature of the LSPRs. On the other hand, the microscale nature of the superlattice requires a very large simulation domain.¹⁰¹ To tackle this challenge, two approaches are generally taken: (i) an effective medium theory approach which neglects the nanoscale effects and periodicity to focus on the overall optical properties of the supercrystal, and (ii) the use of a unit cell with periodic boundary conditions which neglects the overall habit of
the supercrystal to focus on nanoscale behaviors. While the former fails to account for local effects such as interparticle coupling and anisotropy effects as shown by Mirkin and coworkers,³⁶, ⁸⁵ the latter approach, which is used for the calculation of the photonic band structures, fails to describe the photonic properties rising from finite-size and microscopic effects.³⁶, ¹⁰²

Here, we theoretically investigate the photonic properties of 3D periodic microscale arrangements of plasmonic nanoparticles that form supercrystals, also known as plasmonic supercrystals, 3D superlattices, or plasmonic photonic crystals. As I have discussed in Chapter 2, finite-size effects can be unraveled using non-traditional computational methods. In this work, we proposed a theoretical and computational approach for calculating the photonic band structures of finite supercrystals, investigated the effect of the habit (i.e., crystal shape) and size on its photonic properties, and conducted a detailed comparison with standard electrodynamic simulations of the optical properties (reflectance) of these finite superlattices. We also conducted a direct comparison with photonic band structure calculations of truly infinite superlattices. Although we applied this new approach to 3D Au NPs supercrystals, this method is far more general and can straightforwardly be applied to finite 1D and 2D periodic structures of many kinds of building block materials.

### C. Methods

In order to calculate the optical properties of the supercrystals, we use FDTD package from Ansys-Lumerical.⁵⁶ We have considered multiple cases: Infinite and finite crystals, different crystal habits (for the finite case), lattice parameters, and building blocks (nanospheres and nanorods) to determine the photonic band structure and the reflectance of each case. Different methods were used depending on the system and are described below. For all these cases, the cubic unit cell (cf. Figure 10a, real lattice) has been chosen because it is the simplest lattice to analyze. However, it is important to note right away that the proposed method also applies to other unit cell symmetries. In this work, we used an optical excitation in the UV-vis-NIR region (0.83–3.3 eV in energy or 375–1500 nm in wavelength), and which includes the LSPR frequency of the Au NPs (~2.3 eV or 550 nm for Au NPs of 80 nm in size). For calculating the band diagram, only the first-order Brillouin zone is considered. We defined 80 different wavevectors within the first Brillouin zone of the cubic reciprocal lattice (Figure 10a) running over the high symmetry path described by the high symmetry points: $\Gamma$ (center of the Brillouin zone), $X$ (center of a face), $M$ (center of an edge), and $R$ (corner point); thus, requiring one simulation per wave vector, $k$. The crystal building blocks are taken as either Au NPs of 80 nm in diameter or Au nanorods with a size of 40×80 nm immersed in an aqueous medium. Gold is modeled with a frequency dependent refractive index taken from Johnson and Christy's experimentally tabulated data,⁵⁸ and the aqueous medium with a constant refractive index $n_m = 1.33$. Two lattice parameters, $a$, were considered: 100 and 130 nm, corresponding to interparticle distances of 20 and 50 nm, respectively. These lengths, NPs sizes and interparticle distances, were selected in order to maintain the computational resources relatively low without compromising the main goal of the present contribution which is to demonstrate the new approach to calculate the photonic band structure.

#### 1. Infinite Case

Infinite cubic superlattices of Au NPs (Figure 10a) were simulated to determine their photonic band structures. The unit cell was defined by Bloch boundary conditions in x-, y-, and z-directions, used to define the periodicity of the system copying the fields at one edge of the simulation region and re-injecting them at the other edge, ensuring the phase correction of
the fields *via* the Bloch conditions (Figure 10b). The simulation domain (i.e., unit cell) was discretized with a uniform mesh of 2.5×2.5×2.5 nm. Twenty electric dipole sources were randomly distributed within the unit cell, each one with a random phase and a random spatial orientation (defined by the angles $\theta$ and $\phi$) to excite all possible photonic modes within the system. This approach will ensure that all modes are excited even if one dipole is located at a node. We also randomly distribute twenty electric field monitors within the cell to record the time evolution of the electromagnetic field. Using multiple time monitors at different locations ensures all modes are captured. Convergence was assured *via* the auto-shutoff parameter in the Ansys-Lumerical software, which was set at $10^{-7}$ to allow for the electromagnetic fields to propagate for 500 fs. Finally, the photonic band structure was determined from the photonic density of states by applying the fast Fourier transformation (FFT) to the collected electric field in the time domain for each k-vector, where the intensity refers to the squared amplitude of the electric field. The optical signal was previously filtered by a Gaussian function centered at 250 fs and width of 65 fs to ignore initial and final effects.

## 2. Finite Case

Many of the parameters described for the infinite case were used for the present one and are, therefore, omitted here. Three supercrystal structures of different habits were simulated to obtain their photonic band structure and reflectance spectra: Cube, rhombic dodecahedron, and sphere; all of them were structured with Au NPs and ordered in a simple cubic lattice arrangement (Figure 10c). The dimensions of each supercrystal, $L$, were 2 and 3.5 µm (for $a = 100$ and 130 nm, respectively) in length for the cube, face-to-face length for rhombic dodecahedron, and diameter for the sphere. The simulation domain was delimited with standard perfectly matched layers (PMLs, 12 layers) to absorb the electromagnetic fields at the boundaries in the x-, y-, and z-directions. The simulation domain was taken with a size equal to the length of the supercrystal, $L$, plus the maximum source wavelength, i.e., 3.5 µm (for $a = 100$ nm) and 5 µm (for $a = 130$ nm) and was discretized with a non-uniform mesh. A refined mesh of 4×4×4 nm was applied to the Au NPs. The key physical parameters of each system are summarized in Table 2. For the photonic band structure calculations, 6 electric dipole sources were randomly distributed and oriented within a unit cell and distributed to the other unit cells by applying the Bloch conditions (Figure 10b) to the dipole phase in order to couple all these sources within the crystal lattice depending on their spatial positions (Figure 10d). A thousand electric field monitors were also randomly distributed within the entire supercrystals. The photonic band structures were also determined following the same process as for the infinite case. For the reflectance calculations, the supercrystal was excited using an optical plane wave source (total-field scattered-field, TFSF source) with normal incidence at one end and a planar monitor was located at the other end of the supercrystal to record the reflected signal during the simulation.

<table><caption>**Table 2.** Key physical parameters for each system: supercrystal overall size ($L$), unit cell volume ($V_{cell}$), nanoparticle volume ($V_{NP}$), and filling fraction ($F = V_{NP}/V_{cell}$).</caption><thead><tr><th>Supercrystals</th><th><i>L</i><br>[µm]</th><th><i>V</i><sub>cell</sub><br>[nm<sup>3</sup>]</th><th><i>V</i><sub>NP</sub><br>[nm<sup>3</sup>]</th><th><i>F</i></th></tr></thead><tbody><tr><td>Spherical NPs, <i>a</i> = 100 nm</td><td>2.0</td><td>10<sup>6</sup></td><td>2.68×10<sup>5</sup></td><td>0.268</td></tr><tr><td>Spherical NPs, <i>a</i> = 130 nm</td><td>3.5</td><td>2.20×10<sup>6</sup></td><td>2.68×10<sup>5</sup></td><td>0.122</td></tr><tr><td>Nanorods, <i>a</i> = 100 nm</td><td>2.0</td><td>10<sup>6</sup></td><td>8.38×10<sup>4</sup></td><td>0.084</td></tr></tbody></table>
### 3. Effective Medium Model

Maxwell-Garnett theory<sup>87</sup> was used to calculate the effective refractive index ($N_{eff} = \epsilon_{eff}^{1/2}$) of a simple cubic superlattice ($a = 100$ and $130$ nm) composed of Au NPs (80 nm in diameter) in an aqueous medium. Briefly, the effective dielectric permittivity, was defined as $\epsilon_{eff} = \epsilon_m(2\alpha F + 1)/(1 - \alpha F)$, where $\alpha = (\epsilon_{Au} - \epsilon_m)/(\epsilon_{Au} + 2\epsilon_m)$ and $F = V_{NP}/V_{cell}$ is the filling factor (i.e., volume fraction). $\epsilon_{Au}$ and $\epsilon_m = n_m^2$ are the dielectric permittivities of Au and the surrounding medium, respectively, and $V_{NP}$ and $V_{cell}$ are the volumes of the Au NPs and the unit cell, respectively. The resulting $N_{eff}$ was used to determine the dispersion relation, $\omega = kc/N_{eff}$, of light in this medium, and for a slab of 2 µm of thickness using FDTD in order to compare them with that obtained from the finite supercrystal method.

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 10.** (a) Real space simple cubic lattice (left) formed with Au NPs showing the dipole sources (blue arrows) and field monitors (red crosses) within the unit cell used for the infinite case simulation and the first Brillouin zone of the reciprocal cubic lattice (right) showing the high symmetric path (Γ–X–M–R–Γ) used for the band structure calculations. (b) Block condition. (c) Supercrystals with cubic, rhombic dodecahedral, and spherical habits, used in the finite case simulations made up of Au NPs. (d) Simulation layout of the finite supercrystal (9 unit cells are shown) consisting in Au NPs (yellow spheres), electric dipoles (blue arrows), and field monitors (red crosses).</figcaption></figure>

### D. Results and Discussion

First, we calculated the photonic band structure of an infinite cubic lattice for two different lattice parameters, $a$, of 100 nm (interparticle distance of 20 nm, **Figure 11a**) and 130 nm
(interparticle distance of 50 nm, **Figure 11b**) in order to compare them with the finite cases (different habits). We can see that the thin band coming from zero energy at the **Γ** point along the symmetry line **Δ** (between **Γ** and **X**), which is the light dispersion line in the lattice, is bent at ~2.3 eV (i.e., 550 nm) forming a wide band, in the case of *a* = 100 nm. This frequency corresponds to the LSPR of the 80 nm Au NPs, which is confirmed by a simple Mie theory calculation.¹⁰³ Due to the proximity of the NPs, plasmon hybridization occurs, thus producing that broad and intense band in the photonic band structure.¹⁰⁴ When the interparticle distance increases (*a* = 130 nm, **Figure 11b**), the coupling between NPs rapidly decreases, and the light is almost completely absorbed by the individual NPs, resulting in a fading of the photonic band observed in the former band structure (**Figure 11a**). In both cases, there is a bandgap at higher energies due to the light absorption mediated by electronic interband transitions in gold. Using an effective medium approximation (Maxwell-Garnett theory), we have calculated the relation dispersion of the Au NP lattice (*a* = 100 nm) to better visualize the coupling effect of the LSPR and the light.

<figure><img src="image_25.png" /><figcaption>**Figure 11.** Photonic band structures of an infinite cubic lattice composed of 80 nm Au NPs and with a lattice parameter of (a) 100 nm and (b) 130 nm.</figcaption></figure>

Now that the photonic band structure of an infinite Au NPs lattice has been introduced, let us turn our focus on finite-size supercrystals with cubic, rhombic dodecahedral, and spherical habits (**Figure 4**). In general, it is noted that, in all the cases, the photonic band structures seem to qualitatively resemble that of the infinite supercrystal: *i.e.*, exhibiting an intense photonic band around 2 eV due to the light-plasmon coupling and the formation of the collective, hybridized plasmon response. However, we note that new features arise as a result of the supercrystal finite size, and which clearly depend on the supercrystal habit. The first thing to notice is a significant broadening of the various photonic bands calculated for the finite-size supercrystals. Because we are now working outside of the bandgap, the characteristics of a finite periodic structure are a non-monotonous function of the number of periods, which results in surface (*i.e.*, interface) scattering. These losses, in turn, translate into a decrease of the photonic mode lifetime (*i.e.*, increased damping). While these effects have not been investigated for plasmonic-based superlattices, there have been studied in dielectric photonic crystals.¹⁰⁵, ¹⁰⁶

Centering now our attention on the cubic supercrystal with *a* = 100 nm (**Figure 4a**), it is appreciated that multiple dispersion lines are generated at low energies along the symmetry line **Δ** and replicated along the symmetry line **Z** (between **X** and **M**). These new features correspond to Fabry-Pérot photonic modes (FPM), and which occurred from the interference of the reflected waves between the parallel surfaces forming two opposite faces of the cubic supercrystal.³³ In order to confirm the formation of FPM, we calculated the reflectance of the
cubic supercrystal for normal incidence (red arrow in **Figure 5a** inset), i.e., $\mathbf{k} = \mathbf{k}_X$ ($X$ point in the photonic band diagram). The reflectance spectrum, which is shown in **Figure 5a**, exhibits well defined FPMs patterns, previously observed in backscattering measurements in similar nanoparticle assemblies.³⁶ These FPMs are also plotted in the dispersion diagram (**Figure 5b**, magenta crosses) to clarify their connection with the aforementioned bands along the symmetry line $\Delta$. It is worth noting that the FPMs are also coupled to the LSPR (bending and almost complete suppression around the plasmon band). Thanks to this coupling, a bandgap is formed between the two branches approximately at the Au NPs LSPR, termed polaritonic bandgap (PBG).³⁶ The upper branch is quenched due to the electronic interband transitions of gold. This PBG is not formed in the case of an infinite lattice because of the lack of FPM (cf. **Figure 5b**). In the vis-NIR region, the skin-depth of metals such as gold and silver is typically below 25 nm. We also observed that the incident electric field penetrates into the cubic supercrystal up ~ 100–150 nm from the face, where the wave is impinging at the frequency of the upper branch (490 nm). This penetration is due to the absorption from the interband transitions, highly absorbed at the PBG (539 nm) and reflected at the lower branch (650 nm). The same behavior is observed when the supercrystal is rotated by 90° with respect to the incident electric field, i.e., when $\mathbf{k}$ is at the high symmetric point $M$. This behavior is in excellent agreement with previous theoretical work by Ross *et al*.⁸⁶

<figure><img src="image_26.png" /><figcaption>**Figure 12.** Photonic band structures of the finite cubic-lattice supercrystals with (a,b) cubic, (c,d) rhombic dodecahedral, and (e,f) spherical habits. The supercrystals are composed of 80 nm Au NPs with lattice parameters of 100 nm (left panels) and 130 nm (right panels).</figcaption></figure>
Going back to the photonic band structures of the cubic supercrystal (Figure 12a,b), we can see that by increasing the interparticle distance from 20 to 50 nm (Figure 12b), the LSPR remains present but the FPMs are no longer coupled to it; thus leading to the disappearance of the PBG. It can also be noted how the PBG is formed around 550 nm in the band structure of the finite crystal, in contrast with the continuous band in the infinite case. Additionally, we calculated the dispersion relation for a slab with the effective refractive index $N_{\text{eff}}$ of a Au NP supercrystal with a lattice parameter $a = 100$ nm and in aqueous medium. This calculation allows us for verifying the formation and evolution of the FPM bands, which are the only modes present in the slab configuration. This dispersion relation is to be compared with the photonic band structures shown in Figure 12 and which contain a richer variety of photonic features. For this, the FPMs found in the reflectance spectra calculated for the three supercrystals (Figure 14a,c,e) were verified by taking their frequency positions ($\omega = 2\pi c/\lambda$) and applying the resonance condition of a Fabry-Pérot cavity to find the corresponding wavenumber, $k = \omega\text{Re}[N_{\text{eff}}(\omega)]/c$.³⁶ These FPMs fit with the photonic features as showed with the magenta crosses in the dispersion relations (Figure 14b,d,f). Interestingly, at lower energies (i.e., longer wavelengths), the FPM positions calculated with the effective medium theory deviate slightly from those calculated with the full model (Figure 14). This deviation, which has already been observed by Park et al.,³⁶ is attributed to the decreased size parameter, $\pi L N_{\text{eff}}(\lambda)/\lambda$, and is due to the low refractive index and small size-to-wavelength ratio at longer wavelengths. This further illustrates the limitations of the effective medium approach for a quantitative description of finite-size effects such as Fabry-Pérot resonances.

<figure><img src="image_27.png" /><figcaption>**Figure 13.** Mode splitting to plasmon mode energy ratios ($\Omega_R/\omega_0$) at the X point for the finite supercrystals with (a,b) cubic and (c,d) rhombic dodecahedral habits $\hbar\omega_0 \sim 2.3$ eV.</figcaption></figure>
We have also calculated the mode splitting energy ($\hbar\Omega_R$, **Figure 13**) defined as the energy difference between the upper and lower polaritonic bands, in order to quantify the FPM LSPR coupling. We determined the mode splitting to be $\hbar\Omega_R = 0.55$ eV. Although the upper polaritonic band is almost vanished when $a = 130$ nm it is still possible to determine the mode splitting energy to be $\hbar\Omega_R = 0.55$ eV. Considering the LSPR energy to be $\omega_0 = 2.3$ eV, the corresponding normalized coupling strength is $\eta = \Omega_R/\omega_0 \approx 0.24$. While these coupling strengths are much smaller than these obtained in the ultra-strong coupling regime achieved by Mueller *et al.* for Au NP,$^{79}$ it is important to also notice that the gap-to-diameter ratio they have is 0.02–0.12 as opposed to 0.25–0.63 in our work).

<figure><img src="image_28.png" /></figure>

**Figure 14.** Reflectance at the X point ($\mathbf{k} = \mathbf{k}_X$, left panels) with the FPMs indicated by red markers (the red arrow in the inset indicates the incident wavevector). Dispersion relation (right panels) of Au NPs supercrystals with a lattice parameter of 100 nm. (a,b) Cubic, (c,d) rhombic dodecahedral, and (e,f) spherical habits. In the dispersion relations are included the light dispersion in water (cyan solid line), dispersion relation calculated with FDTD of a slab with effective refractive index of the supercrystal (white dashed line), and dispersion relation of the light in the effective medium (blue dashed line) with the FPMs obtained from the reflectance spectrum and applying the resonant condition of a Fabry-Pérot cavity, to find the k value, $k_N = \omega_n \text{Re}[N_{\text{eff}}(\omega)]/c$ (magenta crosses).

The photonic band structure of a rhombic dodecahedral supercrystal (**Figure 12c,d**) is slightly different from that of a cubic supercrystal. While this supercrystal habit appears more
complex, it is nonetheless a system which can be experimentally fabricated.³⁶, ⁷⁸ It is,
therefore, a very relevant morphology to look at. Similarly to the cubic supercrystal, it also
exhibits a PBG, as well as lower and upper polaritonic branches due to the FPM LSPR
coupling, with a splitting energy $\hbar\Omega_R = 0.55$ eV, which is comparable to the values reported
by Park *et al*.³⁶ This coupling almost disappears when the lattice parameter is increased from
100 to 130 nm; however, the splitting energy could be determined to be $\hbar\Omega_R = 0.53$ eV.
These splitting energies are similar to those obtained for the cubic supercrystals because the
filling fractions of both supercrystals are the same ($F = 0.268$). However, the main difference
arises from the emergence of FPM bands along the other symmetry lines, and not only along
$\Delta$ and $Z$ as in the cubic supercrystal. The other directions which sustain these FPM are the
symmetry lines $T$ (between $M$ and $R$) and $\Lambda$ (between $R$ and $\Gamma$). These characteristic FPM
bands are also apparent since the rhombic dodecahedron has parallel surfaces not only along
the $X$-direction, but also along the diagonals of the lattice ($R$ direction; i.e., $\mathbf{k} = \mathbf{k}_R$) causing
Fabry-Pérot interferences. However, it is worth noting that the contrast between these FPMs
and the plasmon band is much smaller than in the cubic supercrystal. The rhombic
dodecahedral supercrystal does not have faces along $\mathbf{k} = \mathbf{k}_M$ and, therefore, it does not
exhibit Fabry-Pérot interferences at the $M$ point. The FPM bands are less intense in the
dispersion diagram of the rhombic dodecahedral supercrystal (Figure 14d) than in the cubic
supercrystal because the Fabry-Pérot interferences occur three times along the same
symmetry line as opposed to only once in rhombic dodecahedron.

Finally, the photonic band structure (Figure 12e,f) and reflectance spectrum (Figure 14e) for the case of a spherical supercrystal presents oscillating patterns similar to those observed in the cubic and rhombic dodecahedral supercrystals. However, since the sphere does not have any parallel flat faces which could sustain FPM, these oscillations are associated with whispering gallery modes.¹⁰⁷, ¹⁰⁸ One can also notice that the spectral features observed in the sphere are narrower than those observed for the FPM in the cube and rhombic dodecahedron, which is another characteristic of WGMs. When considering a spherical optical resonator, such as our sphere supercrystal, characterized by its diameter, $L$, and effective refractive index $N_{\text{eff}}$, the expected spectral positions of the WGM resonances can be found analytically;¹⁰⁸ however, this is beyond the scope of the present study. As it can be seen from the band structure, these bands are present along the entire symmetry path ($\Gamma$-X-M-R-Γ) indicating that these WGM occur independently of the light direction (isotropic), which should come as no surprise considering the spherical symmetry of the supercrystal. Interestingly, the PBG does not appear as clearly in the photonic band structure compared to the two previous cases (Figure 14f), although the reflectance shows the coupling between the WGMs and the LSPR (decrease of the WGMs intensity near the plasmon resonance). While the reason remains unclear, could be related to the small number of Au NPs from the supercrystal surface that are resonating and coupling with the WGMs.

Lastly, we calculated the reflectance spectra (Figure 15a) and photonic band structure (Figure 15b) of a cubic supercrystal composed of isotropically-oriented nanorods with an aspect ratio of 2 (80 nm long and 40 nm wide) to highlight the versatility of the method. These nanorods are also arranged in a simple cubic lattice with lattice parameter $a = 100$ nm, where the long axis of the nanorods is oriented along the symmetry line $\Delta$ (between $\Gamma$ and $X$), as shown in the inset of Figure 15b. Supercrystals composed of elongated nanoparticles have previously been shown to display interesting optical anisotropy/birefringence and polarization dependent properties.³¹, ⁸⁵, ⁸⁶ The reflectance spectra calculated both two different polarizations (Figure 15a), clearly illustrate the polarization dependence of the optical properties of these nanorod supercrystals. Because of the difference in effective refractive index along the two directions, the Fabry-Pérot resonances along these directions
are at different wavelengths. From the band structure (Figure 15), we can clearly see the formation of the PBG between 2.25 eV (longitudinal LSPR) and 2.35 eV (transverse LSPR). More interestingly, we now observe the formation of a new low-dispersion photonic band at lower energy (i.e., longer wavelengths) along the symmetry lines **Δ** and **Z** (i.e., **Γ**–**X**–**M**) that returns to the lower polaritonic branch along the symmetry lines **T** and **Λ** (i.e., **M**–**R**–**Γ**). The appearance of this photonic branch, which dispersion strongly depends on the incident **k** vector clearly illustrates the optical anisotropy of such supercrystal, which has previously been studied by Ross et al.⁸⁶ The strong optical anisotropy of the Au nanorods directly impacts the coupling of the FPMs with the longitudinal and transverse LSPRs, thus leading to coupling/decoupling phases along the various symmetry paths of the photonic band diagram.¹⁰¹ Another clear effect of this anisotropic coupling is the presence of FPM features along the symmetry line **Δ**.

<figure><img src="image_29.png" /><figcaption>**Figure 15.** (a) Reflectance spectra calculated for an incident polarization along the nanorods long axis (blue) and short axis (red) of a finite supercrystal with cubic habit composed of Au nanorods with a length of 80 nm and a width of 40 nm, arranged in a simple cubic lattice with a lattice parameter of 100 nm. (b) Photonic band structure of the same finite cubic supercrystal.</figcaption></figure>

E. Conclusions

Plasmonic-based 3D supercrystals with well-defined habits and with lattice parameters in the deep sub-wavelength regime behave as plasmonic cavity and have, therefore, a strong potential for studies and applications in cavity quantum electrodynamics (QED). While these structures have been shown to have significant optical ability and structural flexibility, the full understanding of their complex optical properties is conditioned by the availability of multiscale computational tools. In this project, we have developed a computational approach, based on the FDTD method that is capable of simultaneously describing nanoscale effects, interparticle coupling, periodicity, and microscale photonic properties associated with the supercrystal morphology. Our approach allows for accurately calculating the photonic band structures from finite-size 3D superlattices of various habits (cubic, spherical, rhombic dodecahedral) and building blocks (spheres and rods). We showed that they differ from the band structures computed for infinite structure and encompass phenomena such as the Fabry-Pérot interferences, whispering gallery modes, the decrease in photonic mode lifetime, and the formation of polaritonic bandgaps.

It is important to notice that, although this approach was applied to 3D supercrystals, it is far more general and versatile. Not only it applies to various morphology, composition, and size
of the building blocks (e.g., cubes, prisms ...), but also to different unit cell symmetries (e.g., fcc, bcc ...), supercrystal habits (e.g., ellipsoidal, pyramidal ...), and dimensionalities of the periodic structures (e.g., 1D chains, 2D arrays).

I anticipate that this new approach will open the door to the study of other types of finite-size periodic structures, such as multi-composition periodic superlattices for the study of exciton photon coupling,¹⁰² and novel magneto-optical properties,³⁵ which can pave the way to new directions in QED and quantum optics.⁷³ Finally, this approach opens up new opportunities for the understanding of fundamental properties of periodic nanostructures as it does not only compare to optical spectroscopy results but also to EELS as shown by García de Abajo and Blanco,³² and as discussed in Chapter 2.
4

# Raman Energy Density in the Context of Acousto-plasmonics

## A. Related Publications

1. Large, N.; Saviot, L.; Margueritat, J.; Gonzalo, J.; Afonso, C. N.; Arbouet, A.; Langot, P.; Mlayah, A.; Aizpurua, J. Acousto-plasmonic hot spots in metallic nano-objects. *Nano Lett.* **2009**, *9*, 3732–3738.

2. Large, N.; Huntzinger, J.-R.; Aizpurua, J.; Jusserand, B.; Mlayah, A. Raman-Brillouin electronic density in short-period superlattices. *Phys. Rev. B* **2010**, *82*, 75310–75318.

3. Montaño-Priede, J. L.; Mlayah, A.; Large, N. Raman energy density in the context of acoustoplasmonics. *Phys. Rev. B* **2022**, *106*, 165425.

## B. Background

As discussed in the previous chapter, there has recently been a great deal of effort made to improve the performance of plasmonic-based applications. Another strategy has been to couple LSPs to other types of elementary excitations, including photonic modes,¹⁰⁹ magnetic modes,¹¹⁰ and excitons.³⁸ To date, this strategy has been very successful as it led to the emergence of new fundamental concepts and technological applications. Most of the studies focusing on metallic nanoparticles consider them as static bodies at rest. However, LSPs can be temporally modulated by acoustic vibrations, also known as acoustic phonons, naturally present in the nanomaterial and surrounding environment. Very recently, there has been a renewed interest on the modulation effect of the LSPs by (opto)mechanical modes and elastic waves.⁹¹, ¹¹¹⁻¹¹⁶ Such high-frequency (GHz-THz) modulation can be used to increase detection and sensing accuracy of nano-objects and molecules. For instance, these vibrating nanoparticles, acting as simple optomechanical nanoresonators, can be seen as nanoscale analogs to quartz crystal microbalances. Another work by Yoo *et al*. showed that epsilon near zero (ENZ) cavities can be used to achieve ultra-strong coupling between ENZ gap plasmon modes and optical phonons.¹¹⁷

Acoustic-phonons Raman scattering is a very effective, high-precision, and non-invasive technique for nanometrology. It uses the acoustic vibrations as local probes to determine sizes and distances at the nanoscale with very high accuracy. However, the effectiveness of this technique is directly related to the knowledge of the interaction mechanisms between electrons (plasmons) and phonons (acoustic vibrations). Therefore, it is necessary to know what are the vibrational modes sustained by a given nanostructure (amplitude, frequency, symmetry) and know how they interact with the LSPs to completely determine the Raman efficiencies and selection rules.¹¹⁸
Here, we investigate the dynamic properties of metallic nanoparticles by focusing on the interaction between confined acoustic vibrations and localized surface plasmons.⁹¹, ¹¹⁹ The dynamic properties of LSPs are responsible for the transient optical absorption modulation, that can be observed in time-resolved transient absorption experiments,¹¹¹, ¹¹², ¹¹⁶ and for the acoustic SERS effect.¹¹⁵ The study of the acousto-plasmonic properties gives unprecedented insight into the fundamental interactions between elementary excitations at the nanoscale and opens up new ways to probe the mechanical properties of nanostructures using optical spectroscopy. Particularly, the acousto-plasmonic interaction in metallic nanoparticles has been studied in low-frequency Raman scattering experiments, in the 2 to 50 cm⁻¹ spectral region, where light is scattered by acoustic vibrations in the nanostructure.¹¹⁵ Acousto-plasmonic driven Raman scattering can be described in three steps: i) incident photon absorption, ii) phonon emission/absorption by the lattice, and iii) scattered photon emission.⁹¹, ¹²⁰ The photon phonon interaction is not direct, but rather occurs via the acousto-plasmonic coupling, i.e., photon electron phonon interaction, which results in the direct modulation of the LSP by the acoustic vibration.¹²¹ When the metal nanostructure is optically excited close to its LSP resonance (LSPR), it leads to a strong acousto-plasmonic interaction. Similarly to the concept of near electric field hot spots in plasmonic nanostructures, the sites where acoustic vibrations produce large modulations of the plasmon induced localized electric field at the nanostructure surface are called acousto-plasmonic hot spots.⁹¹

Here, we described the resonant Raman scattering process using a new approach based on a single direct acousto-plasmonic interaction process described by Fermi's golden rule.¹¹⁹ Within this framework, we introduce a new physical quantity, namely, the Raman energy density (RED). Similarly to the Raman-Brillouin electronic density (RBED) introduced in semiconducting nanostructures,¹⁰⁰, ¹²² the RED allows for studying and monitoring the acousto-plasmonic Raman scattering in the near-field. We modeled the acousto-plasmonic interaction by implementing vibrational dynamics (resonant ultrasound method, RUS),¹²³ into electrodynamic calculations (discrete dipole approximation, DDA).¹²⁴ We used this methodology to compute the acoustic Raman spectra and investigate the interaction between LSPs and (an)isotropic acoustic vibrations, thus leading to the full determination of Raman efficiencies and selection rules. We show that the RED, which can be mapped in the near-field for each vibration mode, correlates the far-field Raman scattering to the local acousto-plasmonic hot spots, provides a deeper understanding of the Raman scattering process (e.g., Raman selection rules), and serves to study the acousto-plasmonic coupling of complex nanostructures.

## C. Theoretical Framework

### 1. Raman Energy Density (RED)

Fermi's golden rule gives the rate at which transitions take place between an initial and a final state of the system consisting of interacting plasmons and phonons in our case. In this framework, the Raman scattering process is described as a single step transition between plasmon and phonon states represented by the occupation numbers $|n_{pl}, n_{ph}\rangle$ ($n_{pl}$ and $n_{ph}$ being the occupation of the plasmon and phonon states). The probability of transition per unit time associated with the Raman scattering process as

$$ \mathcal{R}_{fi}^{(1)} = \frac{2\pi}{\hbar} |\langle n_{pl,f}, n_{ph,f} | \hat{H}_{\text{int}} | n_{pl,i}, n_{ph,i} \rangle|^2 \delta(E_f - E_i), \quad (9) $$

where $E_{i(f)}$ is the energy of the initial (final) plasmon-phonon state. For phonon emission (Stokes Raman scattering) $E_f - E_i = \hbar\omega_{pl,f} + \hbar\omega_{ph} - \hbar\omega_{pl,i}$, where $\omega_{pl}$ and $\omega_{ph}$ are the
plasmon and phonon angular frequencies. $\hat{H}_{\text{int}}$ is the Hamiltonian for the acousto-plasmonic interaction and is expressed as

$$ \hat{H}_{\text{int}} = \sum_{k} \int_{V} -E_{s,k}^{*} \frac{\partial P_{i,k}}{\partial k} u_{m,k}^{*} dV \hat{a}_{s}^{\dagger} \hat{a}_{i} \hat{b}_{m}^{\dagger}, \quad (10) $$

where $\hat{a}^{\dagger}$ (respectively $\hat{a}$) and $b^{\dagger}$ (respectively $\hat{b}$) are the plasmon and phonon creation (respectively annihilation) operators and $E_{s,k}$, $P_{i,k}$, and $u_{m,k}$ are the k-component of the classical electric field, polarization, and displacement vector of phonon mode m, respectively. In previous works,¹⁰⁰, ¹²¹, ¹²² the Raman scattering process was described using a third-order quantum perturbation theory approach. In this three-step process, the surface plasmon is excited by the incoming photon and decays into another surface plasmon state via emission or absorption of confined vibration, and finally a scattered photon is emitted by this plasmon state. However, it is important to notice that our approach uses Fermi's golden rule to describe a one-step transition between two excited plasmon states because the plasmonic fields calculated through electrodynamic simulations already account for the optical excitation and emission processes.

In the absence of any vibration (i.e., nanoparticle at rest), no transition can occur since the LSP states form a set of orthogonal eigenstates. When the nanoparticle vibrates, the LSP polarization vector is modulated, thus enabling dipolar transitions between LSP states. The interaction matrix element between the confined acoustic vibrations and the LSP states is expressed for phonon emission as¹¹⁹, ¹²¹

$$ \langle 1_f 1_m | H_{\text{int}}^0 | 1_i 0 \rangle = - \int \mathbf{E}_f(\mathbf{r}) \cdot \delta_m \mathbf{P}_i(\mathbf{r}) dV, \quad (11) $$

where $\mathbf{P}_i(\mathbf{r})$ is the LSP induced polarization excited in the Raman scattering process (initial LSP state) and modulated by the confined acoustic vibration mode m [$\delta_m \mathbf{P}_i(\mathbf{r})$] and $\mathbf{E}_f(\mathbf{r})$ is the local electric field associated with the final LSP state, which gives rise to the scattered light experimentally detected. Using this interaction matrix element and defining the Raman energy density (RED) as $U_{R,m}(\mathbf{r}, \omega_i, \mathbf{k}_i, \omega_f, \mathbf{k}_f) = -\mathbf{E}_f(\mathbf{r}) \cdot \delta_m \mathbf{P}_i(\mathbf{r})$, we can express the transition rate [Eq. (9)] as

$$ \mathcal{R}_{fi}^{(1)} = \frac{2\pi}{\hbar} \left| \int U_{R,m}(\mathbf{r}, \omega_i, \mathbf{k}_i, \omega_f, \mathbf{k}_f) dV \right|^2 \delta(E_f - E_i). \quad (12) $$

The RED is a complex local energy density that gives rise to the Raman scattering. It contains the excitation, interaction, and emission steps of the Raman scattering process and is in J/m³. It is important to note that there are two fundamental differences between the RED and the RBED previously introduced for semiconductors.¹⁰⁰, ¹²² First, the RED combines electric fields and polarization modulation into a product homogeneous to an energy density while the RBED considers electronic wave functions to give an effective electronic density. Second, contrary to its semiconductor analog, the RED includes the vibrational component through $\delta_m \mathbf{P}_i(\mathbf{r})$. The RED can be divided into two terms:

$$ U_{R,m}(\mathbf{r}, \omega_i, \mathbf{k}_i, \omega_f, \mathbf{k}_f) = -\epsilon_0 \delta_m \chi(\mathbf{r}) \mathbf{E}_f(\mathbf{r}) \cdot \mathbf{E}_i(\mathbf{r}) - \epsilon_0 \chi(\mathbf{r}) \mathbf{E}_f(\mathbf{r}) \cdot \delta_m \mathbf{E}_i(\mathbf{r}), \quad (13) $$

where $\mathbf{E}_i(\mathbf{r})$ and $\mathbf{E}_f(\mathbf{r})$ are the classical local electric fields associated to the LSP before and after the emission of the acoustic vibrations, respectively. $\chi$ is the electric susceptibility of the metallic nanoparticle and is here defined using the Drude model

$$ \chi(\omega_i, \mathbf{r}) = \chi^{\text{ib}}(\omega_i) - \frac{\omega_p^2}{\omega_i^2 + i\omega_i\gamma(\omega_i, R)}, \quad (14) $$
where $\chi$ is dependent on the nanoparticle radius, $R$, through the size corrected Drude damping $\gamma(\omega_i, R) = \gamma_0 + g_s(\omega_i)v_F/R$, where $v_F$ is the Fermi velocity, $\gamma_0$ is the Drude damping in the bulk material, and $g_s$ is obtained from a quantum treatment of the electron surface scattering and is only weakly dependent on $\omega_i$.

### Volume Mechanism

The first term in Eq. (13) describes the modulation of the electric susceptibility by the acoustic vibrations through $\delta_m\chi(\mathbf{r})$ and is called the deformation potential coupling mechanism or volume mechanism.¹²⁰ A change in the nanoparticle volume results in a change in the electronic band structure via the deformation potential and hence in the electric susceptibility. This particular mechanism has been the subject of a recent study by Saison-Francioso et al., who focused on the computational study of shape, electron density, and interband transition effects through the deformation potential mechanism.¹¹³ The modulation of the intraband transitions and electric susceptibility is negligible in the visible range (where the LSP is located) because their excitation lies in the infrared range. Whereas, the interband transition modulation contributes to the volume mechanism as they occur in the UV-visible range, close to the LSP. The volume mechanism (first) term in Eq. (13) reads as

$$U_{R,m}^{\text{VM}}(\mathbf{r}) = -\epsilon_0\chi(\mathbf{r}) \left[ \frac{V_{\text{DP}}}{\hbar\omega_i - \hbar\omega_{\text{ib}}} \nabla \cdot \mathbf{u}_m(\mathbf{r}) \right] \mathbf{E}_f(\mathbf{r}) \cdot \mathbf{E}_i(\mathbf{r}), \quad (15)$$

where $\hbar\omega_{ib}$ is the interband transition threshold, $V_{DP}$ is the deformation potential, and $\nabla \cdot \mathbf{u}_m(\mathbf{r})$ is the divergence of the displacement field $\mathbf{u}_m(\mathbf{r})$ associated with the acoustic vibration mode $m$ and obtained from the elastodynamic calculations.

### Surface Mechanism

The second term in Eq. (13) corresponds to the surface orientation coupling mechanism. Contrary to the volume mechanism,¹²⁰ the surface mechanism is the dominant contribution to the Raman scattering process in metals because of the strong plasmonic field localized at the nanoparticle surface.¹¹⁹, ¹²¹ As it is well known, the LSPs strongly depend on the nanoparticle shape. Therefore, the polarization induced by the LSP at the nanoparticle surface experiences a modulation by the acoustic vibrations through a change in the nanoparticle shape. A simple approximation to account for this surface orientation mechanism consists in a geometrical framework making use of the geometric factors in the nanoparticle polarizability tensor. Here, however, we rigorously describe the surface mechanism [second term in Eq. (13)] in terms of the difference between the polarization fields of the final and initial LSP states.

## 2. Raman Scattering Spectrum

The Stokes Raman scattering spectrum of a metallic nanoparticle, derived from Fermi's golden rule for the acousto-plasmonic coupling [Eq. (12)], is expressed as

$$I_R = \sum_m \frac{2\pi}{\hbar} \left| \int U_{R,m}(\mathbf{r}) dV \right|^2 \frac{\Gamma_m}{(\hbar\omega_f + \hbar\omega_m - \hbar\omega_i)^2 + (\Gamma_m/2)^2} \left[ \frac{1}{e^{(\hbar\omega_m/k_B T)} - 1} + 1 \right] \quad (16)$$

where the summation runs for all the acoustic vibration modes. $\Gamma_m$ and $\omega_m$ are the spectral linewidth and eigenfrequencies, respectively, of the acoustic vibration mode $m$. The $\delta(\hbar\omega_f + \hbar\omega_m - \hbar\omega_i)$ function is defined using a Lorentz distribution to account for the homogeneous broadening. The term $[e^{(\hbar\omega_m/k_B T)} - 1]^{-1}$ is the Bose-Einstein population factor considering a Stokes process. As can be seen in Eq. (16), integrating the local Raman energy density yields the far-field low-frequency acoustic Raman scattering spectrum from the plasmonic nanoparticle.
<table><caption>**Table 3.** Gold parameters used for the RED and Raman spectrum calculations.</caption><thead><tr><th colspan="3">Electrodynamics Simulation (DDA)</th><th colspan="3">Elastodynamics Simulation (RUS)</th></tr><tr><th>Parameter</th><th>Value</th><th>Reference</th><th>Parameter</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>χ<sup>ib</sup></td><td>J&C database</td><td>58</td><td>ρ</td><td>19,700 kg/m<sup>3</sup></td><td>125</td></tr><tr><td>ħω<sub>p</sub></td><td>9.01 eV</td><td>126</td><td>v<sub>L</sub></td><td>3240 m/s</td><td>125</td></tr><tr><td>γ<sub>0</sub></td><td>0.07 eV</td><td>126</td><td>v<sub>T</sub></td><td>1200 m/s</td><td>125</td></tr><tr><td>g<sub>s</sub>v<sub>F</sub></td><td>0.915 eV nm</td><td>126</td><td>ν</td><td>0.42</td><td>125</td></tr><tr><td>V<sub>DP</sub></td><td>0.8 eV</td><td>127</td><td></td><td></td><td></td></tr><tr><td>ħω<sub>ib</sub></td><td>2.4 eV</td><td>126</td><td></td><td></td><td></td></tr></tbody></table>

## D. Model and Methods

### 1. Computational Model

To demonstrate and illustrate the concept of RED, we calculated the RED and computed the Raman spectrum of a free-standing Au spherical NP of radius $R = 2.5$ nm (Figure 16a), immersed in water (refractive index $n = 1.33$), and optically exited at its LSPR ($\lambda = 520$ nm, Figure 16b). The incident optical excitation is a plane wave of unitary amplitude, polarized in the y-direction, and traveling along the x-direction.

In our model, we neglect any mechanical interactions between the NP and the surrounding aqueous medium. However, it is important to note that the mechanical coupling between a more rigid environment and the nanoparticle may need to be taken into consideration. Furthermore, the heat generated by the NP under optical excitation has a negligible effect on the vibration frequencies and the LSPR, and will be neglected in our calculations.

### 2. Elastodynamics: Resonant Ultrasound Spectroscopy (RUS) Method

RUS is a method based on the measurement of the frequencies of the free vibrations of a solid object in order to determine the elastic tensor of the material of which the object is made. It has been successfully applied to systems ranging from small metallic clusters to macroscopic objects.123 The vibrations are calculated using Rayleigh-Ritz variational approach. The displacement vector field associated with each acoustic vibration mode is calculated assuming an elastic continuous medium and using the xyz algorithm proposed by Visscher et al. to solve the Navier-Stokes hydrodynamic equations.123, 128 RUS uses the longitudinal and transverse sound velocities in the material and expands the displacement vectors onto a basis $Φ_ζ = x^i y^j z^k$, where $ζ = (i, j, k)$ is the function label. Here, we used $0 ≤ i + j + k ≤ 20$ in order to achieve a good numerical convergence.123 The Rayleigh-Ritz approach takes the dynamic problem into a generalized eigenvalue problem where all the eigenmodes obtained are orthonormalized.123 The values used for the Au density ($ρ$), longitudinal and transverse sound velocities ($v_L$ and $v_T$, respectively), and Poisson's ratio ($ν$) are quoted in Table 3.

### 3. Electrodynamics: Discrete Dipole Approximation (DDA)

Once the vibration modes have been calculated in RUS, we use DDA (DDSCAT v7.3 package)124 to calculate the local electric field at the LSPR wavelength for the particle at rest (initial state, Figure 16b) and for the nanoparticle deformed by a given acoustic vibration mode. In all the DDA simulations, we used an interdipole distance $d = 0.03$ nm to discretize the nanoparticle and achieve numerical convergence; the total number of dipoles within the
nanoparticle was around $3\times10^6$. The physical parameters used for gold in Eq. (14) are provided in **Table 3**.

<figure><img src="image_30.png" /><figcaption>**Figure 16.** Calculated absorption spectrum of a 5 nm AuNP in water showing the dipolar LSPR at 520 nm. Inset: 3D color map and (001) cross-sectional plane of the ENF at the LSPR.</figcaption></figure>

## E. Results and Discussion

### 1. Acoustic Vibrations

To compute the Raman scattering spectrum of a AuNP undergoing acoustic vibrations and the RED associated with each vibration, we first calculate the surface displacements of a free-standing, homogeneous, isotropic, and elastic Au nanosphere induced by the acoustic vibration modes. In Lamb's original paper, the acoustic vibration modes were classified into torsional and spheroidal modes.¹²⁹ Torsional modes do not induce any change in the material density, which implies that the divergence of the displacement is $\nabla \cdot \mathbf{u}_m = 0$, thus leading to $U_{R,m}^{\text{VM}} = 0$. Moreover, because of the absence of shape change induced by such modes, we also have $U_{R,m}^{\text{SM}} = 0$. Therefore, because they do not contribute to either the surface or the volume mechanism, these torsional modes are Raman inactive and will thus not be considered in this work. On the other hand, spheroidal modes induce changes in the NP shape and/or volume. These modes are labeled $S_{\ell m}^n$, where the integers $n$, $\ell$, and $m$ denote the harmonic ($n = 1$ being the fundamental), the angular momentum number, and its z-component azimuthal number, respectively. The vibrational density of states is discrete and the mode eigenfrequencies are given by $\omega_{\ell,n}[\text{cm}^{-1}] = \xi_{\ell,n}v_L/2R$, where $\xi_{\ell,n}$ is a mode-dependent coefficient. **Figure 17** shows the surface displacements associated with the fundamental breathing mode ($S_{00}^1$; **Figure 17b**) and the fivefold degenerated fundamental quadrupole modes ($S_{2m}^1$ with $m = 0, \pm 1, \pm 2$; **Figure 17c-g**).

Because gold exhibits a strong elastic anisotropy, we have also calculated the surface displacements associated with the six equivalent breathing and quadrupolar vibration modes from an anisotropic Au nanosphere. The irreducible representations of the anisotropic vibration modes have been determined based on the $O_h$ point group character table.¹³⁰ As a result of the elastic anisotropy, the spheroidal breathing modes transform into an $A_{1g}$ vibration (**Figure 17h**). Similarly, the anisotropy induces a partial degeneracy lift of the fivefold spheroidal quadrupolar mode which splits into two $E_g$ (**Figure 17i-j**) and three $T_{2g}$ (**Figure 17k-m**) degenerated vibrations; the latter can be projected onto spheroidal Lamb modes.¹³⁰
<figure><img src="image_31.png" /><figcaption>**Figure 17.** (a) AuNP at rest. Top: NP deformed by the isotropic breathing $S_{00}^1$ ($d_0$) mode (b) and the fivefold degenerated quadrupole $S_{2m}^1$ ($d_1 - d_5$) modes (c-g). Bottom: NP deformed by the anisotropic breathing $A_{1g}$ mode (h), the quadrupole $E_g$ [$E_{g(a)}$ and $E_{g(b)}$] modes (i-j), and the quadrupole $T_{2g}$ [$T_{2g(a)}$, $T_{2g(b)}$, and $T_{2g(c)}$] modes (k-m).</figcaption></figure>

The eigenfrequencies and the labels we use for each vibration mode are presented in **Table 4**. Because of the weak mechanical coupling between the nanoparticle and the surrounding aqueous environment, the nanoparticle can be considered as free standing. Therefore, we also compare their eigenfrequencies obtained using RUS to analytical results from Lamb theory.¹²⁹

<table><caption>**Table 4.** Eigenfrequencies of the isotropic and anisotropic vibration modes calculated using RUS and compared to analytic results from Lamb theory.</caption><thead><tr><th colspan="2">Lamb theory</th><th colspan="2">Isotropic (RUS)</th><th colspan="2">Anisotropic (RUS)</th></tr><tr><th>Mode</th><th>Frequency<br>$S_{lm}^n$ (cm<sup>-1</sup>)</th><th>Mode</th><th>Frequency<br>(cm<sup>-1</sup>)</th><th>Mode</th><th>Frequency<br>(cm<sup>-1</sup>)</th></tr></thead><tbody><tr><td>$S_{00}^1$</td><td>20.8</td><td>$d_0$ (Figure 17b)</td><td>20.8</td><td>$A_{1g}$ (Figure 17h)</td><td>20.7</td></tr><tr><td>$S_{20}^1$</td><td>7.1</td><td>$d_1$ (Figure 17c)</td><td>7.1</td><td>$E_{g(a)}$ (Figure 17i)</td><td>5.0</td></tr><tr><td>$S_{2\pm1}^1$</td><td>7.1</td><td>$d_2$ (Figure 17d)</td><td>7.1</td><td>$E_{g(b)}$ (Figure 17j)</td><td>5.0</td></tr><tr><td>$S_{2\mp1}^1$</td><td>7.1</td><td>$d_3$ (Figure 17e)</td><td>7.1</td><td>$T_{2g(a)}$ (Figure 17k)</td><td>8.0</td></tr><tr><td>$S_{2\pm2}^1$</td><td>7.1</td><td>$d_4$ (Figure 17f)</td><td>7.1</td><td>$T_{2g(b)}$ (Figure 17l)</td><td>8.0</td></tr><tr><td>$S_{2\mp2}^1$</td><td>7.1</td><td>$d_5$ (Figure 17g)</td><td>7.1</td><td>$T_{2g(c)}$ (Figure 17m)</td><td>8.0</td></tr></tbody></table>

It is clear from **Figure 17** that the directions of maximum deformation are mode dependent. The breathing mode ($d_0$) is isotropically deformed in all directions due to its pure radial nature, thus preserving the spherical shape throughout the vibration period (Figure 17b). On the other hand, the anisotropic breathing mode ($A_{1g}$) has maximum deformations in the [100] directions, i.e., along the x-, y-, and z-directions, thus breaking the spherical symmetry (Figure 17h). For the sake of presentation, we have chosen a cross-sectional plane of the NP which contains the direction of maximum deformation for each vibration mode. These planes of maximum deformation are also chosen to contain the polarization direction (y-direction) of
the optical excitation used in the electrodynamic calculations. Although more pronounced for the breathing modes, the difference on the surface deformations between isotropic and anisotropic vibration modes is clearly visible.

## 2. Plasmonic Near-Field

In order to compute the RED, we have calculated the electric near-field (ENF) of the AuNP at rest excited at its LSPR (Figure 16) with a plane wave propagating in the x-direction and polarized along the y-direction (i.e., initial state). The inset in Figure 16 shows the ENF spatial distribution at the surface of the NP and in the (001) plane (xy-plane); it clearly shows the dipolar nature of the LSP and a maximum electric field enhancement $|E/E_0| \approx 5$.

<figure><img src="image_32.png" /></figure>

**Figure 18.** ENF induced by the dipolar LSP and modulated by the isotropic (top) and anisotropic (bottom) vibration modes. The cross-sections are taken where the maximum deformation is located; the (xyz) planes are indicated in each panel. (a) NP at rest, (b) isotropic breathing $S_{00}^1$ ($d_0$) mode, (c-g) fivefold degenerated isotropic quadrupoles $S_{20}^1$ ($d_1 - d_5$), (h) anisotropic breathing $A_{1g}$ mode, (i,j) anisotropic quadrupole $E_g$ [$E_{g(a)}$ and $E_{g(b)}$] modes, and (k-m) anisotropic quadrupole $T_{2g}$ [$T_{2g(a)}, T_{2g(b)},$ and $T_{2g(c)}$] modes.

**Figure 18** shows the NEF spatial distributions on the surface of the NPs deformed by the vibration modes as well as in the cross-sectional planes of maximum deformation. As discussed earlier, this figure shows how the ENF is distinctively and spatially modulated by each individual vibration mode. However, it is important to notice that the amplitude of the modulation depends on the mode. This, will ultimately reflect on the Raman efficiency of the mode as we will discuss further.

The ENF distribution of the AuNP deformed by the isotropic breathing mode (Figure 18b) does not appear to change in comparison with that of the NP at rest (Figure 18a). This absence of modulation from the isotropic breathing mode is due to (i) the conservation of the spherical shape and (ii) the size change smaller than the optical wavelength and negligible retardation effects. Consequently, according to Eq. (13), the surface mechanism contribution to the RED will be negligible. On the other hand, the anisotropic breathing mode (Figure
**18h)** appears to have the largest ENF enhancement between the analyzed modes due to the strongest localization of surfaces charges at the NP poles (along the y-direction). The same occurs for the *d*₁ and *E*<sub>g(a)</sub> vibration modes (**Figure 18c,i**), where the deformation of the NP is also along the applied electric field, yielding a little more enhancement compared to the NP at rest. For the other modes, the ENF enhancement is imperceptibly modulated because the surface displacements of the vibrating NP occur in directions other than that of the applied electric field polarization.

**3. Raman Energy Density (RED)**

We have calculated the partial RED associated with the volume and surface mechanisms as well as the total RED [Eq. (13)] (in neV/nm³) for each isotropic and anisotropic vibration mode. It should be noted that the spatial region between the surface of the nanoparticle at rest and the deformed surface was not considered and set to zero due to the lack of physical realism. Indeed, the surface location cannot be defined with a precision better than the interdipole spacing in the DDA simulations.

<figure><img src="image_33.png" /><figcaption>**Figure 19.** Real part of the surface $U_{R,m}^{SM}$ (a,d), volume $U_{R,m}^{VM}$ (b,e), surface+volume $U_R$ (c,f) components of the RED associated with the isotropic (a-c) and the anisotropic (d-f) breathing modes.</figcaption></figure>

**Figure 19** displays the spatial distributions of the real part of $U_{R,m}^{SM}$, $U_{R,m}^{VM}$, and $U_{R,m}$ for the isotropic (**Figure 19a-c**) and anisotropic (**Figure 19d-f**) breathing modes. These maps show the spatial extent of the coupling between the acoustic vibrations and the dipolar LSP yielding inelastic light scattering in and around the NP. They represent the acousto-plasmonic local source of the Raman scattering measured in the far-field. The partial RED associated with the deformation of the NP surface ($U_{R,m}^{SM}$) by the $d_0$ mode is localized only on the NP surface; there is no acousto-plasmonic coupling inside the NP. When anisotropy is included ($A_{1g}$ mode), $U_{R,m}^{SM}$ is enhanced thanks to the localization of the LSP at the deformed surface and to the acousto-plasmonic coupling occurring inside the NP. On the other hand, $U_{R,m}^{VM}$ shows a radially increasing acousto-plasmonic coupling from the surface to the center of the NP, which is mainly induced by the larger expansion (or compression) inside the NP than underneath the surface.<sup>111</sup> The volume contribution to the RED, as seen in panels b and e, arises from the local variations of the electric susceptibility induced by the vibration modes through the deformation potential mechanism [Eq. (15)]. This result is in good agreement with the local volume variations of the dielectric permittivity produced by the NP vibration modes and which have been recently predicted by Saison-Francioso *et al.*<sup>113</sup>
The total RED ($U_R$, **Figure 19c,f**) clearly appears to be more spatially localized in the case of the $A_{1g}$ mode compared to the $d_0$ mode. These RED hot spots correspond to regions where the NP simultaneously exhibits a strong surface deformation (from the vibration) and a strong induced ENF (from the LSP). Nevertheless, it is important to notice that the Raman activity of each acoustic vibration mode is obtained by spatially integrating the RED [Eq. (16)]. As a result of this spatial integration, positive and negative quantities may cancel each other out, thus leading to a decrease or to a complete cancellation of the Raman response, despite the presence of RED hot spots. Based on this observation, the $d_0$ mode is expected to lead to a larger Raman efficiency than the $A_{1g}$ mode, because the former has no negative contribution to the RED.

<figure><img src="image_34.png" /><figcaption>**Figure 20.** Real part of the total RED, $U_{R,m}$, associated with the isotropic vibration modes [$d_0$ (a), $d_1$ (b), $d_2$ (c), $d_3$ (d), $d_4$ (e), $d_5$ (f)] and anisotropic vibration modes [$A_{1g}$ (g), $E_{g(a)}$ (h), $E_{g(b)}$ (i), $T_{2g(a)}$ (j), $T_{2g(b)}$ (k), $T_{2g(c)}$ (l)].</figcaption></figure>

The same analysis applies to the other isotropic and anisotropic modes. **Figure 20** shows the total RED, $U_{R,m}$, calculated for each isotropic and anisotropic mode. For most of the fivefold degenerated quadrupolar modes, the contribution of $U_{R,m}^{SM}$ in the total RED is larger than $U_{R,m}^{VM}$ in both the isotropic and anisotropic cases. Nevertheless, as discussed above, because of cancellation of the negative and positive contributions to the RED, one expects a lower Raman efficiency for $A_{1g}$. As for example, it could be expected that the contribution of the $d_3$, $d_4$, and $d_5$ modes (**Figure 20d-f**) to the Raman scattering will be low, also for $T_{2g}$ modes (**Figure 20j-l**). On the other hand, it is also expected that the contribution of the $d_1$ and $E_{g(a)}$ modes will be large owing to their high RED intensity inside the NP, mostly from the $U_{R,m}^{SM}$ contribution. Although $d_2$ and $E_{g(b)}$ modes have positive values (**Figure 20c,i**, respectively), their contribution to the Raman scattering is likely to be low due to their weak RED values.
4. Acoustic Raman Spectra

**Figure 21a-c** displays the Raman spectra calculated for each isotropic vibration mode and considering the surface and volume mechanisms as well as their sum, respectively.

The *d*₁ mode gives the strongest Raman response because its maximum displacement is parallel to the incident electric field polarization. On the other hand, *d*₀ and *d*₂ have very weak contributions through the surface mechanism because the former leads only to small changes of the NP shape and the latter exhibits the maximum deformation in the direction perpendicular to the incident polarization. The only mode that leads to a strong contribution to the Raman efficiency through the volume mechanism is *d*₀. When the NP undergoes a *d*₀ mode, it experiences a global expansion (or compression) unlike the other modes, which on the contrary experience simultaneous compression and expansion at different locations on the NP surface, and thus resulting in the cancellation of the Raman efficiency.

<figure><img src="image_35.png" /><figcaption>**Figure 21.** Raman spectra calculated using Eq. (16) for the surface (a,d), volume (b,e), and both (c,f) mechanisms of the isotropic (a-c) and anisotropic (d-f) vibration modes.</figcaption></figure>

On the other hand, the Raman spectra calculated for each anisotropic vibration mode and considering the surface and volume mechanisms as well as the sum of the two mechanisms are presented in **Figure 21d-f**, respectively. Similarly to the isotropic case, the quadrupolar modes produce a strong Raman intensity through the surface mechanism while the breathing mode contributes through the volume mechanism. However, a few differences can be noticed: the T₂g modes Raman peaks are shifted to higher frequency and are less intense than their isotropic counterparts *d*₃, *d*₄, and *d*₆ due to the anisotropy in sound velocities.

As a consequence, the Raman spectrum of an isotropic AuNP has principally two peaks, which correspond to the acousto-plasmonic coupling of the fivefold degenerated quadrupolar vibration modes *via* surface mechanism and of the breathing mode *via* the volume mechanism (**Figure 22**, blue). The Raman spectrum of an anisotropic AuNP has three peaks.
The two lower-frequency peaks are due to the quadrupolar $E_g$ and $T_{2g}$ modes coupled to the LSP *via* the surface mechanism while the higher-frequency peak is due to the breathing $A_{1g}$ mode which is Raman active *via* the volume coupling mechanism (Figure 22, red). As previously discussed, the positive and negative contributions of the RED associated with the anisotropic breathing mode cancel each other out, thus leading to a decrease in Raman intensity and resulting in the $d_0$ mode being more intense than the $A_{1g}$ mode. As for the quadrupolar modes, it is important to note that we recover the intensity relationship between the isotropic and anisotropic modes. Indeed, the sum of the intensities of the isotropic quadupolar modes (i.e., $\sum_{k=1}^5 d_k$) is equal to the sum of the intensities of the anisotropic quadrupolar modes (i.e., $\sum_{k=a}^c E_{g(k)} + T_{2g(k)}$), with a relative error of 1.5$\times 10^3$.

<figure><img src="image_36.png" /><figcaption>**Figure 22.** Total Raman spectra calculated using Eq. (16) for the isotropic and anisotropic modes.</figcaption></figure>

In both the isotropic and anisotropic cases, the quadrupolar $d_1$ and $E_{g(a)}$ modes are the most intense as previously discussed. These results are in excellent quantitative agreement with several experiments performed by on small AuNPs.<sup>115, 131</sup> In particular, the measured Raman intensity of the $T_{2g}$ modes is weaker as compared to the $E_g$ modes, in agreement with our theoretical calculations (Figure 22, red). Once again, the difference in intensity between these two anisotropic quadrupolar mode symmetries lies in their RED components. Indeed, as shown in Figure 20, the real parts of the RED associated with the $T_{2g}$ modes (panels j-l) exhibit positive and negative contributions, which tend to cancel out in Eq. (12). On the other hand, the real parts of the RED associated with the $E_g$ modes (panels h-i) are predominantly either positive or a negative.

**5. Raman Selection Rules**

Finally, we compute the Raman spectra using Eq. (16) for various incident and scattered polarization configurations to determine the Raman selection rules. We demonstrate the contribution of each mode in the situation where an analyzer (polarization filter) parallel or perpendicular to the incident electric field is used before the Raman detector (Figure 23). First, we see that the Raman intensities have slightly decreased when a parallel analyzer is used in both the isotropic and anisotropic cases. Furthermore, when using a perpendicular analyzer (cross configuration), the breathing mode contribution becomes negligible and only a few quadrupolar contributions persist. It is well known from experiments that the breathing mode of a spherical nanoparticle is totally polarized, whereas the quadrupolar modes are only
partially polarized. This can easily be understood when looking at the various RED components. Due to its spherical symmetry, the modulation produced by the isotropic breathing mode, through the surface mechanism, has no component perpendicular to the incident polarization. However, it is important to notice that the breathing mode does not completely vanish, due to the weaker, yet contributing, volume component of the RED.

<figure><img src="image_37.png" /><figcaption>**Figure 23.** Raman spectra calculated using Eq. (16) for the parallel (a,c) and cross (b,d) configurations of the isotropic (a,b) and anisotropic (c,d) modes.</figcaption></figure>

Conversely, the quadrupolar modes induce a RED with components along the incident field direction. Their individual symmetries, which is given by relative orientation of the major deformation axis with respect to the incident polarization, dictates their relative Raman intensities (Figure 23). These observations are in excellent agreement with the theoretical work by Duval,¹¹⁸ who showed that for the Raman transition that excites a spherical breathing mode, the polarizations of incident and scattered fields are parallel only, while they can be parallel or perpendicular in the case of the quadrupolar modes. Our theoretical approach also explains the selections rules observed experimentally by Tripathy et al. in the framework of acoustic SERS,¹¹⁵ Interestingly, we can note that the experiments conducted by Tripathy et al. revealed a large number of harmonics for the various modes. While we limited our work to the fundamental modes (n = 1) only, our approach can easily include harmonics to describe the experimental observations. For the harmonics (n > 1), two important factors will impact the Raman intensity. First, the inner deformation associated with the presence of displacement nodes within the volume will directly impact the volume component of the RED. Second, the lower surface displacement amplitudes will result in weaker modulation of the ENF at the nanoparticle surface, thus directly impacting the surface component of the RED. Both of these effects lead to a decrease of the Raman efficiency with increasing overtones, making them challenging to observe.¹²⁰ While an atomistic description of the nanoparticle vibration eigenmodes is traditionally required,¹³² our RED framework provides a new route to fully describe such high-order harmonics in the acoustic Raman spectra.
## F. Conclusions

We introduced the Raman energy density (RED) as a new theoretical tool for the interpretation of resonant Raman scattering mediated by LSPs in metallic nanostructures. Analogous to the LDOS, it represents the electromagnetic energy density, which is excited by the Raman probe and modulated by the acoustic vibrations of the nanostructure. It is a local quantity that can be mapped in the near-field region to provide a clear picture of the acousto-plasmonic interaction, which gives rise to the inelastic light scattering measurable in the far-field. It allowed us to calculate for the first time the Raman selection rules taking into account the surrounding medium and the modulation of the polarization in this medium. We showed that the RED can accurately predict the Raman active modes and the ones forbidden by the Raman selection rules for both isotropic and anisotropic modes; in excellent quantitative agreement with other theoretical frameworks<sup>118, 132</sup> and experiments.<sup>115, 132</sup> The RED provides significant insight about the origin of the plasmon-vibration coupling (volume vs. surface mechanisms).

It is worth noting that the RED framework goes well beyond plasmonic nanoparticles and can be applied to all types of nanomaterials including semiconductors and dielectrics. For instance, dielectric nanoparticles sustain Mie resonances, which induce local electric fields. In such case, similar electrodynamic simulations can be performed and coupled to the elastodynamic calculations to compute the resulting acoustic Raman scattering spectra and associated RED. Semiconductor quantum dots (QD) can also be described using this approach. While the local electric field induced by excitons is much weaker, leading to a small surface contribution, the volume deformation is expected to strongly modulate the QD electronic density of state, thus resulting in a strong volume contribution. Another aspect that can also be integrated into this framework concerns the mechanical properties of the surrounding medium, which affect the nanoparticles vibrational properties.<sup>133</sup> Our approach can capture and provide insight into the mechanical and electromagnetic coupling through the RED and computed Raman spectra.

The RED constitutes a novel and unique tool for the exploration of physical effects such as acoustic SERS,<sup>115</sup> vibrational transfer,<sup>116</sup> vibrational hybridization,<sup>134</sup> nonlinear mixing,<sup>135</sup> and other optomechanical phenomena between neighboring NPs, substrates, and surrounding media. It constitutes a step towards the realization of nanophononic platforms capable of probing simultaneously electronic, thermal, and phononic states.
# Plasmonic Photothermal Heating Effects in Nanostructures

## A. Related Publications

1. Bustamante, G.; Carrizales, K.; DeLuna, F.; Large, N.; Ye, J. Y. Fabrication and characterization of thermo-responsive gold nanorod assemblies. In *Proc. SPIE*, **2018**, 10507, pp 105070B-1.

2. Lopez, P.; Mayer, K., Large, N, Synthesis, characterization, and computational modeling of polyelectrolyte-coated plasmonic gold nanorods for photothermal heating studies. In *Proc. SPIE*, **2019**, 11082, pp 110822W-1.

3. Manrique-Bedoya, S.; Abdul-Moqueet, M.; Lopez, P.; Gray, T.; Disiena, M.; Locker, A.; Kwee, S.; Tang, L.; Hood, L.; Feng, Y.; Large, N.; Mayer, K. M. Multiphysics modeling of plasmonic photothermal heating effects in gold nanoparticles and nanoparticle arrays. *J. Phys. Chem. C* **2020**, *124*, 17172-17182.

4. Wang, R.; Li, J.; Rigor, J.; Large, N.; El-Khoury, P. Z.; Rogachev, A. Yu.; Kurouski, D. Direct experimental evidence of hot-carrier-driven chemical processes in tip-enhanced Raman spectroscopy (TERS). *J. Phys. Chem. C* **2020**, *124*, 2238-2244.

5. García-Rosas, C. M.; Medina, L. A.; Lopez, P.; Large, N.; Reyes-Coronado, A. Magneto-plasmonic biocompatible nanorice. *J. Nanopart. Res.* **2021**, *23*, 144.

6. Mondal, S.; Montaño-Priede, J. L.; Tu Nguyen, V.; Park, S.; Choi, J.; Hoang Minh Doan, V.; Mai Thien Vo, T.; Hung Vo, T.; Large, N.; Kim, C. S.; Oh, J. Computational analysis of drug-free silver triangular nanoprism theranostic probe plasmonic behavior for in-situ tumor imaging and photothermal therapy. *J. Adv. Res.* **2022**, *41*, 23-38.

7. Rigor, J.; Kurouski, D.; Large, N. Plasmonic heating effects in tip-enhanced Raman spectroscopy (TERS). *J. Phys. Chem. C* **2022**, *126*, 13986-13993.

## B. Background

Light-nanoparticle interactions, as well as light-tissue and nanoparticle-tissue interactions, can be studied using computational modeling to further understand the effect of light diffusion in tissue and the temperature change induced by the heating of nanoparticles when excited with light. The window opened by computational modeling enables the optimization of therapies like plasmonic photothermal therapy (PPTT) as the size and shape of AuNPs can be easily modified in search of the type of particles that will enhance the thermal effects within a target. Other computational models, such as the one presented by Baffou *et al.* implement shape factors to simplify the calculation of the power absorbed by non-spherical
nanoparticles and their subsequent heat generation.136, 137 Another method, proposed by Hogan et al.,138 studies the effect of multiple scattering in heat generation. This is due to the use of nanoshells with a relatively large size (~150 nm), which possess a scattering dominated extinction cross-section.

Here, we took advantage of the tools provided by computational modeling and compare the effects of size and shape by creating computational models of different types of nanoparticles, more specifically, gold nanorods, gold bipyramids (GBPs) and gold nanospheres (GNSs). The light nanoparticle interaction can be explored using these computational models, as the power absorbed by the particles, when excited by an incident light wave, can be calculated from the solutions of Maxwell's equations. Furthermore, the subsequent nanoparticle heating due to the absorption of light and its effect in temperature change of the surrounding medium (water) can be assessed by using the power absorbed as a heat source to solve the heat equation. Our computational model uses a novel approach to calculating power absorbed by the nanoparticles and their subsequent heat generation, as it incorporates a full electrodynamics simulation to solve the light nanoparticle interaction and couple the results into a heat transfer model without implementing shape simplifications. Additionally, we focused on nanoparticles whose extinction cross-section is dominated by absorption and exhibit relatively little scattering, making multiple scattering effects less significant in our model. These simulations were compared with analytical solutions obtained with Mie theory and with experimental results performed in the group of Dr. Katie Mayer (UTSA).

### C. Method

Computational models of different nanoparticles were based on transmission electron microscopy (TEM) images of GNRs, and GBPs. Each one of the models was setup so that different sizes and aspect ratios could be analyzed by modifying specific parameters. In the case of the GNS, this parameter is the radius. For the GNR, the parameters include the length and width. Lastly, the modifiable parameters for the GBP are the length, width, and tip radius.

#### 1. Optical Absorbance

We used the finite element method (FEM, Comsol's RF module)139 to solve Maxwell's equations. An absorbing boundary, or perfectly matched layer (PML), with a thickness of $λ/2$ was implemented to truncate the computational domain (Figure 24).

<figure><img src="image_38.png" /><figcaption>**Figure 24.** Computational model used in (a) the electrodynamic simulations, showing the PMLs around a nanorod and (b) the heat transfer simulations, showing a water sphere around a bipyramid.</figcaption></figure>
In addition to the absorbing property of the PMLs, and to reduce the impact of the discretization, a scattering boundary condition was assigned to the outermost surface, meaning that the waves reaching the external boundary of the PMLs would not travel back into the simulation domain.

<figure><img src="image_39.png" /><figcaption>**Figure 25.** Calculated optical efficiencies: extinction (blue), absorption (green), and scattering (red). The numerical solutions (solid lines) are compared to the analytical solutions from Mie theory (markers) for (a) 150 nm Au nanosphere in air and (b) 50 nm Au nanosphere in water.</figcaption></figure>

Illumination at the LSPR of the nanoparticles could result in either light scattering or light absorption. The amount of scattered and absorbed light was determined using the cross-section calculations and summed up to obtain the extinction cross-section of the nanoparticles (Figure 25). However, for the scope of this work, we focused on the light absorption and the resulting localized heating. The absorption cross-section of the particle can be calculated from Poynting's theorem as

$$ \sigma_{\text{abs}} = \frac{W_{\text{abs}}}{S_{\text{in}}}, \qquad (17) $$

where $S_{\text{in}}$ is the Poynting vector magnitude and where the power absorbed by the NP is

$$ W_{\text{abs}} = \int_V Q_h dV. \qquad (18) $$

Here, $Q_h$ is defined as the total power dissipation density in W/m³ (i.e., the total losses of the system) based on the calculated electromagnetic fields. We obtained the total power absorbed by integrating over the volume, $V$, of the GNP.

## 2. Absorbed Power Density

The absorbing properties of nanoparticles play an important role in photothermal heating. At the LSPR wavelength, the light absorbed is maximized and so is the photothermal heating effect. To calculate the photothermal heating generated by a spherical nanoparticle, Eq. (19) yields an accurate approximation based on the absorption cross-section, the volume of the particle, and the irradiance of the light source.

$$ Q_{\text{abs}} = \frac{\sigma_{\text{abs}}I_s}{V} \qquad (19) $$

This approach, however, yields inconsistent results when treating complex geometries such as the bipyramid or the nanorod; in these cases, a more general approach is undertaken as shown in Eq. (20):
$$Q_{\text{abs}} = \frac{W_{\text{abs}}}{V}, \qquad (20)$$

with $Q_{\text{abs}}$ being the power per unit volume generated by the nanoparticle (in W/m³), and considering that the effect of other losses due to atomic effects are negligible. It is clear to see how the power absorbed by the GNP is readily transformed into heat.

### 3. Heat Transfer Calculation

COMSOL's Coefficient Form PDE module was used for our heat transfer model, due to its versatility for adding or neglecting different effects in the simulations. With a heat source, $Q_{\text{abs}}$, and temperature, $T$, the steady-state heat equation reads

$$-\nabla(\kappa\nabla T) = Q_{\text{abs}}, \qquad (21)$$

where $\kappa$ is the diffusion coefficient (in W/m·K). Equation (21) is where the absorption calculations from the electrodynamic model feed into the heat transfer model. Initially, $Q_{\text{abs}}$ is computed using FEM [Eq. (20)]. This value is then used as the heat source in Eq. (21).¹³⁷

For anisotropic nanoparticles such as the GNR, their optical absorption is polarization dependent. The RF module was used to calculate the average heat power across all possible orientations and use this value as the heat generated by a particle in any orientation or even in an array where nanoparticles are randomly oriented.

Considering that we are not working with propagating waves anymore, the truncation method (PML) is no longer required, and so the size of the domain is not restricted. A 3-µm sphere was placed around the particle to act as surrounding medium (water) as shown in Figure 25b.

## D. Results and Discussion

### 1. Single Au Nanoparticles

The computational model was verified by comparing the results obtained computationally with the analytical results obtained following Mie theory. Since the analytical solution is available only for spherical nanoparticles, we used Mie theory to calculate the spectra of various GNSs across several wavelengths, and with different surrounding media. These analytical solutions were then compared with the numerical simulations obtained when using a computational model with the same surrounding media and a spherical geometry of the same dimensions. Figure 25 shows an excellent agreement between the computational and the analytical results for spheres of two sizes: a 150 nm sphere in air and a 50 nm sphere in water. The average relative error is 2.3% and 4.2% for the nanoparticle in air and water, respectively. It can be noticed that, as expected, a smaller particle has a higher absorption cross-section, whereas the scattering cross-section is more predominant in the larger sized particle.

Furthermore, we also compared results to UV-vis absorbance spectra obtained from samples of GNRs and GBPs synthesized in the group of Dr. Mayer. The LSPR wavelengths were found to be similar between the calculated and measured spectra for both nanoparticle types. However, due to the inhomogeneity in the samples (particles of different sizes and shapes, as well as additional byproducts and aggregates resulting from the synthesis process), the breadth of the peaks observed in the experimental spectra were different from the ones observed in the simulations.

Since these AuNPs can be used as plasmonic photothermal therapy agents, we designed them so they exhibit high absorptive properties in the NIR (~ 800 nm). For this purpose, six different computational models of nanoparticles were compared: two rods, two bipyramids, and two spheres. The GNR and GBP samples synthesized in Dr. Mayer's lab exhibit
resonance wavelengths appropriate for NIR implementation. The spherical models were generated based on the volume measured from the TEM of the Au nanorods and the Au bipyramids. The radii of the modeled spheres were chosen so that their volume would be equal to the volume of the modeled GNR and GBP ($R_{eq}$ = 34 nm in the case of the GBP, and $R_{eq}$ = 22 nm in the case of the GNR). It is worth noting that both width and endcap radius remain unchanged for the two nanorod models ($W$ = 26 nm; $R$ = 11 nm) as well as for the two bipyramid models ($W$ = 58 nm; $R$ = 8 nm). For this reason, the dimensions column of **Table 5** reports only the length, $L$, of the GNR and GBP and the diameters of the GNS. Depending on the size and shape of the nanoparticle, it exhibits either the highest absorbance or the highest scattering when the longitudinal modes are excited at the resonance wavelength.

<table><caption>**Table 5.** Nanoparticle dimensions, their LSPR wavelength, and calculated absorbed power densities.</caption><thead><tr><th>Nanoparticle Type</th><th>Dimensions [nm]</th><th>Dipolar LSPR [nm]</th><th>Absorbed power density at 0° [mW/μm³]</th><th>Absorbed power density at 90° [mW/μm³]</th><th>Avg. absorbed power density [mW/μm³]</th></tr></thead><tbody><tr><td>Nanorods</td><td>91</td><td>807</td><td>933.89</td><td>0.33</td><td>311.52</td></tr><tr><td>Nanorods</td><td>94</td><td>820</td><td>777.68</td><td>0.33</td><td>259.44</td></tr><tr><td>Bipyramids</td><td>153</td><td>807</td><td>248.77</td><td>0.45</td><td>83.22</td></tr><tr><td>Bipyramid</td><td>158</td><td>820</td><td>234.46</td><td>0.44</td><td>78.45</td></tr><tr><td>Nanosphere</td><td>68</td><td>536</td><td>0.80</td><td>0.80</td><td>0.80</td></tr><tr><td>Nanosphere</td><td>44</td><td>528</td><td>0.67</td><td>0.67</td><td>0.67</td></tr></tbody></table>

With our model, the maximum light absorption is measured when the wave's propagation vector, **k**, is perpendicular to the longitudinal axis of a non-spherical particle, and the polarization of the wave, **E**, is parallel to the longitudinal axis. The opposite configuration yields the minimum amount of power absorbed: **k** parallel to the longitudinal axis of the GNP and **E** perpendicular. The orientation-dependent nature of the power absorption was accounted for by calculating a mean power absorbed, $Q_{abs}^{avg}$, as the average between three configurations of light propagation: the first configuration is where **E** and the longitudinal axis of the particle are parallel, $Q_{abs}^{0°}$, and the other two configurations, grouped together as they yield the same value, are where **E** and the longitudinal axis are perpendicular, $Q_{abs}^{90°}$. This average absorbed power density can be written as

$$Q_{abs}^{avg} = \frac{Q_{abs}^{0°} + 2Q_{abs}^{90°}}{3}, \qquad (22)$$

and is reported in **Table 5** for each nanoparticle.

Solving Maxwell's equations in the light-nanoparticle system allowed for the calculation of the electric field enhancement induced by the plasmon excitation in the GNPs. The magnitude of both absorbed power and field enhancement also depends on the light intensity used. Based on the typical values found in the literature,¹³⁷ we chose an irradiance of 1 mW/μm² and a wavelength of 808 nm. Given the size difference between the laser beam spot and the nanoparticle, a plane wave model was implemented and the amplitude was calculated using the relationship between the wave's amplitude and irradiance. The effect of size and shape in the electric field enhancement for a single GNR, GBP, and GNS can be observed on the left column of **Figure 26** and a volume map of the total power dissipation
density (power absorbed by the same particles) is shown on the center column of **Figure 26**. This power density was then utilized in Eq. (20) and then ultimately in Eq. (22).

<figure><img src="image_40.png" /><figcaption>**Figure 26.** Left (a-c): Electric field magnitude under a plane wave illumination of 808 nm wavelength and 1 mW/µm² irradiance for (a) a 91 nm nanorod (b) a 153 nm bipyramid, and (c) a 68 nm nanosphere. Center (d-f): Power absorbed by (d) a 91 nm nanorod (e) a 153 nm bipyramid, and (f) a 68 nm nanosphere under the same illumination. Right (g-i): Steady-state temperature maps of (g) a 91 nm GNR, (h) a 153 nm GBP, and (i) a 68 nm GNS. Scale bars are 50 nm.</figcaption></figure>

The results reported in both **Table 5** and **Figure 26** indicate that the nanorods absorb the highest amount of optical power when under illumination of a 808 nm laser. More specifically, the 91 nm long GNR absorbs more power than the 94 nm GNR due to the close proximity between the particle's resonance wavelength and laser's wavelength. Furthermore, it is observed that the bipyramid models show a considerable amount of absorbed power, despite being over 50% larger than the rods. Thus, we may assume that both the rod and the bipyramid geometries could be good candidates for plasmonic heating applications. However, size limitations inherent to the particular application need to be taken into consideration. Additionally, shape effects are observed as the individual spherical models do not absorb a significant amount of power despite having the same volume as the GNR and GBP models; therefore, are not suitable for plasmonic heating in the NIR region due to their low power absorption levels at these wavelengths.

We have shown that both the GNRs and the GBPs absorb significant amounts of power in the NIR region, now we will measure the change in temperature achieved by each particle due to the amount of absorbed power under 808-nm laser illumination. The average power absorbed by the individual particles will be used as a heat source in the heat transfer model [Eq. (21)].
The maximum temperatures reached by the six nanoparticles used in this work, along with the average absorbed power used for the calculations, are reported in **Table 6**. The heat maps of the three nanoparticles exhibiting the highest temperature reached (within each type) are shown in right column of **Figure 26**.

<table><caption>**Table 6.** Average absorbed power density per unit volume for different types of nanoparticles and the maximum temperature reached.</caption><thead><tr><th>Nanoparticle Type</th><th>Dimensions [nm]</th><th>Avg. Absorbed power Density [mW/µm³]</th><th>Max. Temperature [°C]</th></tr></thead><tbody><tr><td>Nanorod</td><td>91</td><td>311.52</td><td>101.1</td></tr><tr><td>Nanorod</td><td>94</td><td>259.44</td><td>91.7</td></tr><tr><td>Bipyramid</td><td>153</td><td>83.22</td><td>80.2</td></tr><tr><td>Bipyramid</td><td>158</td><td>78.45</td><td>78.4</td></tr><tr><td>Nanosphere</td><td>68</td><td>0.80</td><td>37.5</td></tr><tr><td>Nanosphere</td><td>44</td><td>0.67</td><td>37.1</td></tr></tbody></table>

As expected, since the nanorods absorbed the highest amount of power, they also reach the highest temperature out of the three shapes compared, with the 91 nm nanorod being the one that reached the highest temperature (101°C). These results correlate to the trends observed experimentally, where a 655 nm laser was used to heat three samples: GNRs, GBPs, and pure water. After 22 min., the system reached steady state and it was observed that the nanorods induced a temperature change of ~5°C, whereas the GBPs induced a temperature increase of ~2°C. Water alone did not show any temperature increase under laser illumination. Interestingly, both the experiments and the simulations show that the bipyramids can also induce a temperature increase. **Table 6** shows that the smaller GBP reaches a higher temperature than the larger GBP. Similar to the GNR, the smaller bipyramid has a LSPR wavelength closer to the wavelength of the laser, thus the plasmonic effect is enhanced and the subsequent heating is higher. This finding provides a sound argument to the initial thought presented earlier: both GNRs and GBPs can achieve significant heating under NIR illumination, thus making them good candidates for plasmonic heating applications. Because they exhibited the highest temperature increase, we select the 91 nm GNR as the best morphology for plasmonic heating under NIR illumination.

## 2. Periodic Au Nanoparticle Arrangements

Now that we have chosen the 91 nm nanorod as an ideal geometry, we can study ensemble and arrangement effects. Specifically, we designed a 3D lattice of NPs in solution and a 2D array deposited on a substrate. Both the 2D and 3D models measure the effect posed by changing the distance between the particles (i.e., changes in the concentration). We tested two different concentrations: 1 and 22 nM, the former resulted in an interparticle spacing of ~1100 nm, whereas the latter resulted in an interparticle spacing of ~420 nm. These interparticle distances are much larger than the nanoparticle size, so we expect near-field interparticle interactions to be negligible. An approximation of the randomly distributed nanoparticles can be obtained modeling a perfect cubic (3D), or square (2D), lattice array where the distance between particles is controlled by the concentration of the solution. Both arrays are composed of a finite amount of nanoparticles: 27 GNRs in the 3D array and 9 in the 2D array, allowing us to simulate the effects in the micron-scale. Crucially, the heat
power calculated from one nanoparticle in the nanoscale simulation (as described in the previous section) is input into the micron-scale simulation as the new heat source power, thus bridging the two length scales in our computational model. The interparticle spacing in these systems is much smaller than the laser beam spot and the light decay length in the solution, so it can be assumed that each particle absorbs the same amount of power as its immediate neighbors.138 It is worth noting that a randomly distributed array of nanoparticles may enable increased interparticle interactions; future work is underway to evaluate this effect. For this reason, and to account for the random orientation observed experimentally, the average power density absorbed by a single 91 nm long GNR (Table 6) was assigned to each particle in the array and the heat transfer for the different configurations was calculated using a steady state model. Figure 27 shows the temperature maps computed for the 3D model at both concentrations.

<figure><img src="image_41.png" /><figcaption>**Figure 27.** Temperature spatial distribution around an array of 91 nm long nanorods in water with concentrations of 1 nM (top row) and 22 nM (bottom row). (a,c) show the 3D temperature distributions on a quarter section of the arrays. (b,d) show the temperature distributions in the cross-sectional plane of the NPs. The heat sources are set at 312 mW/µm³.</figcaption></figure>

As mentioned above, the average power absorbed by the 91 nm GNR (312 mW/µm³) was assigned to each one of the particles in the array. Each of them contains 27 GNRs arranged on a cubic lattice. The steady state heat transfer simulations from the 3D arrays at different concentrations show that the area where hyperthermia levels are reached (T > 43°C) are similar in size for both systems: the 1 nM system induced hyperthermia levels over a region of 24.4 µm³, whereas the region affected by the heating of the 22 nM system has a volume of 28.7 µm³. Although the volume experiencing the plasmonic heating is similar, the highly concentrated system reached higher temperatures than the diluted one, indicating that a solution with higher concentration of GNPs absorbs greater amounts of light, thus inducing hyperthermia levels much faster while containing the effect within one micron from the center of the array.
For the 2D array, it is important to consider possible interface effects such as nanoparticles on a substrate or at the interface between two different substances.¹³⁶ To check whether interface effects would be significant, one simulation was run with a 2D array of nanorods placed at a water-glass interface and another simulation contained the 2D array surrounded by water only. There was a negligible difference between the optical absorbance between the array at the glass-water interface and the one surrounded by water. For this reason, and to reduce computational time, the heat maps of the 2D arrays at 1 and 22 nM were calculated without substrate and are shown in **Figure 28**. These temperature maps show a behavior comparable to the one observed with the 3D lattice, but the regions where hyperthermia levels are reached are significantly smaller. The 1 nM system induced hyperthermia in the areas surrounding the particles but no further than 400 nm away from each particle. The 22 nM array induced hyperthermic levels within a region of ~ 1.5 µm, smaller than its 3D counterpart but significantly larger than the 2D array at 1 nM. This indicates that the heating effect of GNPs in 2D, similar to the 3D case, is enhanced with increased concentration.

<figure><img src="image_42.png" /><figcaption>**Figure 28.** Temperature distribution for a 2D array of nanorods in water at (a) 1 nM and (b) 22 nM. Each particle was assigned an average heat source power of 312 mW/µm³.</figcaption></figure>

### 3. Random Ag Nanoparticle Ensembles

In collaboration with the group of Dr. Junghwan Oh at Pukyong National University, we invested the potential use of PVA-coated Ag triangular nanoprisms (SNT) in cancer photothermal therapy. To do so we calculated the absorption cross-section of randomly and uniformly distributed and randomly oriented SNT (5 nm thickness), Ag hexagonal nanoplates (NH, 5 nm in thickness), and silver nanospheres (NS) (Figure 29a). The dielectric permittivity tabulated by Palik has been used for silver.¹⁴⁰ We used Gaussian size distributions and relative concentrations based on the experimental size distribution histograms and TEM images. The average size (measured as the radius of their circumcircle), standard deviation, σ, and the number of nanoparticles in the distribution are shown in **Table 7**. The relative concentration, defined as $n_{part}$ is equal to the total count of a particle type divided by the total count of SNT in the histograms, that is, $n_{SNT} = 1$, $n_{NH} = 0.4$, and $n_{NS} = 0.7$.

<table><caption>**Table 7.** Average size, standard deviation, number of nanoparticles, and concentrations used for the electrodynamic simulations of NT, NH, and NS.</caption><thead><tr><th>Nanostructure</th><th>Size [nm]</th><th>σ [nm]</th><th>Number</th><th>n<sub>part</sub>, sim. 1 [NPs/µm<sup>3</sup>]</th><th>n<sub>part</sub>, sim. 2 [NPs/µm<sup>3</sup>]</th></tr></thead><tbody><tr><td>NT</td><td>44.5</td><td>18</td><td>100</td><td>30</td><td>90</td></tr><tr><td>NH</td><td>53.5</td><td>39</td><td>40</td><td>12</td><td>36</td></tr><tr><td>NS</td><td>13.1</td><td>1.9</td><td>70</td><td>21</td><td>63</td></tr></tbody></table>
<figure><img src="image_43.png" /><figcaption>**Figure 29.** (a) Model of a random distribution of silver NT, NH, and NS. (b) Theoretical maximum temperature generated by the SNT (red) and experimental temperature (blue).</figcaption></figure>

We performed electrodynamic simulations to obtain the absorption spectra of each individual nanoparticle ensembles (NT, NH, and NS) with the same concentrations (15 NPs/µm³). The calculated absorption spectra showed that the NT and NH exhibit a wide band in the NIR region while the NS exhibits only a dipolar plasmon mode at 413 nm. Then, two absorption spectra were calculated for a distribution composed of a mixture of the three nanoparticles using the two concentrations specified in **Table 7** ($n_{part}$). At low NP concentration we noted that a broad absorption band appeared in the wavelength range 800–1200 nm due to the plasmonic interaction between the NTs and the NHs. Increasing the concentrations threefold resulted in a blue-shift of the NIR broad band to 700–1100 nm due to the increasing near-field interaction between the nanoparticles. Further increase of the concentration leads to further broadening of the band, similar to the experimental spectrum. Finally, we have calculated the photothermal properties of the SNT ensemble.

**Figure 29b** shows the experimental photothermal steady-state temperatures (blue) at different laser power densities (0.5, 1.0, 1.5, and 2.0 W/cm²) under 808-nm laser irradiation and the theoretical maximum temperatures for the same power densities using an optical excitation at the maximum absorption wavelength. Theoretically, the temperature increases linearly as a function of the applied optical power, as shown in Eq. (19). These results are in agreement with the computational work done by Manrique-Bedoya et al.<sup>141</sup> There is also an excellent quantitative agreement between the temperatures predicted from our computational model and the experimental measured temperatures (within 2% at high power density). The slight deviation from the theoretical model observed at high powers can be explained by the thermal energy lost in the environment (heat dissipation), water volatilization from the cells, and protein denaturation.<sup>142, 143</sup>

4. TERS Tip

Tip-enhanced Raman spectroscopy (TERS) is a spectroscopy technique which possesses single-molecule sensitivity and sub-nanometer spatial resolution. These unique properties are achieved thanks to the extremely high electromagnetic field confinement at the apex of the scanning probe. In collaboration with the group of Dr. Dmitry Kurouski at Texas A&M University, we showed that such strong field confinement can lead to thermal decomposition of the analytes.<sup>44</sup> We demonstrated that the use of an aqueous solvent as tip-sample junction mediator drastically reduces possible molecule degradation. Using a combination of electrodynamic and heat transport simulations, we provide some theoretical insight into the plasmonic heating of the TERS system.
<figure><img src="image_44.png" /></figure>

**Figure 30.** Temperature distribution at the tip-substrate junction and at the surface of the metal when optically excited at 633 nm. (a) Vertical cross-section in the plane (x,z). (b) Horizontal cross-section in the plane (x,y), taken on the substrate (z = 0). The spatial distribution of temperature change, $\Delta T$, is normalized for each of the four cases (Au/air, Au/water, Ag/air, and Ag/water).

**Figure 30** shows the temperature distribution at tip-substrate junction and on the surface of an ITO substrate calculated for an incident optical power of 80 µW, which corresponds to an irradiance $I = 13 \ \mu W/\mu m^2$. First, it can be noted that the whole metallic tip is thermally uniform and is the hottest part of the system, achieving temperature changes as high as $\Delta T_{Ag/air} = 138$ K and $\Delta T_{Au/air} = 69$ K for the Ag and Au tips in air, respectively (**Figure 30a**). When placed in the aqueous solution, the tip temperature is significantly lower than in air; the tip experiences a maximum temperature increase of only a few degrees ($\Delta T_{Au/water} = 2.9$ K, and $\Delta T_{Ag/water} = 5.6$ K). Second, due to its metallic properties, the ITO substrate also heats up slightly, exhibiting a maximum temperature change of the order of $\Delta T_{ITO/air} \sim 2-5.5$ K, with a hot spot immediately underneath the tip location (**Figure 30b**). Similarly to the tip, when the substrate is in contact with water its surface experiences a much lower temperature change. The inset of **Figure 30a** shows a close-up view of the temperature gradient within the 5 nm gap, which is consistent with the theoretical $1/r$ spatial dependence. Finally, it is worth noting that the ratio of maximum temperature change on the tip surface in air and water is $\Delta T_{air}/\Delta T_{water} = 24$; which, unsurprisingly, is the ratio between the thermal conductivities of water and air.

E. Conclusions

We have developed computational models that successfully describe the plasmonic excitation of single gold nanoparticles, the subsequent heat generation due to power absorption from laser illumination and its effects on temperature distribution throughout the surrounding medium. We developed these models in a FEM framework to provide the community with a complete and detailed set of procedures to enable a wide spread of such a type of PPTT modeling. Comparing GNPs of different sizes and shapes confirmed the gold nanorod as the particle shape with the greatest heat power per unit volume, more specifically, the 91 nm GNR was found to have optimal heat generation. Additionally, due to the small size of the nanoparticle and the location of its resonance wavelength in the NIR, the GNR can be used as photothermal agent in applications such as PPTT. Furthermore, the concentrations chosen for
this study represent realistic concentrations of nanoparticles which would be used in medical
applications and it was found that arrays of nanorods can induce hyperthermic levels in the
surrounding medium at low and high concentrations.

We note that our computational model showcases qualitative agreement with the experiments regarding the localized temperature increase both the GNRs and SNTs. Locally concentrated NPs, especially through methods like molecular targeting, can enhance the clustering effects and enable higher temperature increase within the target region as is the purpose of plasmonic heating applications. The computationally-driven theranostic approach using SNTs demonstrated an innovative thought regarding the nanoparticles shape, size, concentration, and composition which could be useful for the prediction of photothermal heat generation in precise nanomedicine applications.

Finally, our simulations show that temperatures, which can reach significant local values of hundreds of kelvins in tip-substrate systems such as in TERS, can cause thermal degradation of analytes present in the probe-sample junction. On the other hand, when TERS is performed in aqueous solutions (L-TERS), the temperature increase of the tip apex is significantly limited to only a few tens of kelvins, ultimately resulting in a more stable TERS signal over time. While these computational models do not account for atomic-scale details (e.g., impurities, crystallinity, grain boundaries) and remain ideal, as compared to experimental systems, they support the use of aqueous solvents as mediators to minimize the thermal degradations of the analytes.
6 |

Perspectives

Moving forward, my work will continue to include a mixture of computational methods, theoretical model development, and close collaborative work with experimental groups for the study of the optical, photothermal, and thermo-optical properties of a wide range of nanomaterials. My future research plans build upon the methods and models I have developed over the past 10 years. Furthermore, all the future research projects, which are presented below, are intertwined and complement each other as I will discuss. For instance, the photothermal model can and will be implemented into the acousto-plasmonic framework. Similarly, the EELS and photonic band structure calculations can and will be compared as they have been shown by García de Abajo to relate to each other. Lastly, although it is not in my short- or mid-term plans, it would be very interesting to look at a periodic system or vibrating nanoparticles and evaluate the effect of the acoustic vibrations on the photonic band structure of such system, in particular if the nanoparticles are mechanically linked through DNA or embedded into a polymer matrix.

A. Acousto-Plasmonics

1. Application of the RED Framework to other Optomechanical effects

With respect to the research area in acousto-plasmonics, the newly-developed framework has the potential to advance, not only the field of plasmonics, but also the field of mechanical properties of nanostructures. Indeed, it opened up new opportunities for the study of a wide range of optomechanical phenomena including: (i) plasmon-mediated vibration transfer, between mechanically uncoupled nanostructures, (ii) beating phenomena in nanoparticle-substrate and coupled nanoparticle systems, and (iii) plasmon-induced optical forces. Understanding the effect of the nanoparticle size, shape, composition, and surrounding environment on the plasmon-phonon coupling can now be done using our new approach. One of the main future objectives, which will be achieved through collaborations, is to directly compare theoretical predictions from the RED formalism to experimental acoustic Raman scattering results.

In interesting perspective lies in a recent work from Dr. Jeremy Baumberg's group at the University of Cambridge, who recently published an experimental study on ultrastrong acousto-plasmonic coupling in nanojunctions composed of a metallic nanoparticle on a metallic substrate.¹¹⁶ In such system, Dr. Baumberg observed the appearance of a highly-localized bouncing acoustic vibration mode. Because of the strong confinement of both electromagnetic field and acoustic vibration in the nanoscale particle-substrate gap, the acousto-plasmonic coupling is extremely strong. While Dr. Baumberg's work uses a simple mechanical model to interpret their experimental results, RED formalism could provide significant insight into the coupling mechanism and the activation of such unusual vibration mode, which I hypothesized to be the result of either vibrational transfer processes between the nanoparticles and the substrate or vibration hybridization.
## 2. Optical Forces Effect on Acoustic Vibrations

I have developed a long-term collaboration with Dr. Adnen Mlayah at the LAAS-CNRS to study these various effects. Dr. Mlayah’s expertise in plasmonics and in acoustic Raman scattering will be useful for experimentally testing the predictions obtained from this new RED theoretical framework. Furthermore, preliminary discussions and preliminary results obtained by my former graduate student (Joel Rigor) during his Chateaubriand Fellowship visit in Toulouse in 2020 lay down the foundation for future studies on plasmon-induced optical forces and their effect on acoustic vibrations (Figure 31).

<figure><img src="image_45.png" /><figcaption>Figure 31. (a) Magnitude of the optical forces (in femto-Newton) for a Ag NP dimer. (b,c) Optical-force-induced deformation of a Ag and Au NP dimer, respectively.</figcaption></figure>

When metallic NPs are optically excited close to their LSPR, optical forces are induced as the result of the presence of strong local electric field induced by the surface plasmons at their surface. These forces contribute to the plasmon-phonon interaction process. Here, I propose to develop a theoretical model to describe and conduct a pioneering study of these effects. Because of the multiphysics and quantum nature of this interaction, we need to develop a unique semi-classical approach that combines electrodynamics, vibrational dynamics, and quantum field theory. This project will be divided into several tasks depicted in Figure 32a. While the work presented in Chapter 3 focused on the effect of the acoustic vibration on the localized surface plasmons, here I will focus on the reciprocal effect, i.e., the effect of the LSPs on the acoustic vibrations. Two of my graduate students (Alex Ferere and John Leblanc) are currently working on these two sides of the plasmon-vibration interaction with a self-consistent method implementation in the Comsol’s FEM framework (Figure 32).

First, we must be able to compute the electric, $E_μ$, and magnetic, $B_μ$, fields operators associated with the NP’s LSP (FEM, Comsol’s RF module). We will then use the canonical quantization of the electromagnetic fields to obtain a rigorous normalization:

$$E_{\mu}(\mathbf{r}) = \sqrt{\frac{\hbar\omega_{\mu}}{2\epsilon_{0}\mathcal{V}_{\mu}}} F_{\mu,E}(\mathbf{r}) \hat{a}e^{-i\omega t} + \text{cst} \quad \text{and} \quad B_{\mu}(\mathbf{r}) = \frac{1}{c}\sqrt{\frac{\hbar\omega_{\mu}}{2\epsilon_{0}\mathcal{V}_{\mu}}} F_{\mu,B}(\mathbf{r}) \hat{a}e^{-i\omega t} + \text{cst}, \quad (23)$$

From these quantized and normalized electromagnetic fields, we can construct the Maxwell’s stress tensor $T_{ij} = T_{ij}^E + T_{ij}^B$, where the electric and magnetic components are:

$$T_{ij}^{E} = \epsilon_{0} \left( E_{i}E_{j} - \frac{1}{2} \delta_{ij} E^{2} \right) \quad \text{and} \quad T_{ij}^{B} = \frac{1}{\mu_{0}} \left( B_{i}B_{j} - \frac{1}{2} \delta_{ij} B^{2} \right). \quad (24)$$