ABSTRACT
Title of dissertation: NUMERICAL STUDY OF PLASMON
RESONANCES IN NANOPARTICLES
Zhenyu Zhang, Doctor of Philosophy, 2007
Dissertation directed by: Professor Isaak D. Mayergoyz
Department of Electrical and Computer Engineering
Surface plasmon resonances in nanoparticles have numerous promising scien-
tiflc and technological applications in such areas as nanophotonics, near-fleld mi-
croscopy, nano-lithography, biosensor, metamaterial and optical data storage. Con-
sequently, the understanding of plasmon resonances in nanoparticles has both fun-
damental and practical signiflcance. In this dissertation, a new numerical technique
to fully characterize the plasmon resonances in three-dimensional nanoparticles is
presented.
In this technique, the problem of determining the plasmon resonant frequen-
cies is framed as an integral equation based eigenvalue problem, and the plasmon
resonant frequencies can be directly found through the solution of this eigenvalue
problem. For this reason, it is computationally more e?cient than other \trial-
and-error" numerical techniques such as the flnite-difierence time-domain (FDTD)
method. This boundary integral equation method leads to fully populated dis-
cretized matrix equations that are computationally expensive to solve, especially
when a large number of particles are involved in the nanostructures. Since the
fully populated matrices are generated by integrals with 1=r-type kernel, this com-
putational problem is appreciably alleviated by using the fast multipole method
(FMM). The boundary integral approach is also extended to the calculation of the
extinction cross sections of nanoparticles, which reveal important information such
as the strength and the full width at half maximum (FWHM) of these resonances.
The numerical implementation of this technique is discussed in detail and numerous
computational results are presented and compared with available theoretical and
experimental data.
Furthermore, metallicnanoshellsforbiosensingapplicationsaswellasnanoparticle-
structured plasmonic waveguides of light are numerically investigated. The integral
equation based numerical technique presented throughout this dissertation can be
instrumental for the design of plasmon resonant nanoparticles and to tailor their
optical properties for various applications.
NUMERICAL STUDY OF PLASMON RESONANCES
IN NANOPARTICLES
by
Zhenyu Zhang
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulflllment
of the requirements for the degree of
Doctor of Philosophy
2007
Advisory Commmittee:
Professor Isaak D. Mayergoyz, Chair/Advisor
Dr. Nail A. Gumerov, Co-Advisor
Professor Ramani Duraiswami
Professor Perinkulam S. Krishnaprasad
Professor Howard M. Milchberg
Professor Thomas E. Murphy
Professor Howard C. Elman
c Copyright by
Zhenyu Zhang
2007
DEDICATION
To My Parents, Jianlin and Xinli.
ii
ACKNOWLEDGMENTS
First and foremost I would like to thank my advisor, Professor Isaak D. May-
ergoyz for introducing me to the fleld of plasmonics and providing me an invaluable
opportunity to work on the challenging and interesting problems over the past four
years. He has always made himself available for help and advice. It has been a great
pleasure to work with and learn from such an extraordinary individual.
I would also like to thank my co-advisors, Dr. Nail A. Gumerov and Professor
Ramani Duraiswami. They have not only taught me the fast multipole method
but also supported me in many ways. Thanks are due to my collaborators Pro-
fessor Donald Fridkin from University of California, San Diego, Professor Gennady
Shvets from the University of Texas at Austin and Professor Giovanni Miano of the
University of Naples Federico II, Italy for numerous enlightening discussions.
I also wish to thank Professor Perinkulam S. Krishnaprasad, Professor Howard
M. Milchberg, Professor Thomas E. Murphy and Professor Howard C. Elman for
agreeing to serve on my dissertation committee and for taking their precious time
reviewing the manuscript.
Furthermore, my friends and colleagues at the University of Maryland have
enriched my graduate life in difierent ways and deserve a special mention: Mr.
Seokjin Kim, Dr. Tao Shen, Dr. Heng-Tung Hsu, Dr. Su Li, Dr. Zhiyun Li, Mr.
Patrick McAovy, Dr. Iulian Nistor, Dr. Mihai Dimian, Dr. Petru Andrei, and Dr.
iii
Chun Tse.
I owe my deepest thanks to my family - my father and mother who have always
been loving and supporting me during this long journey. Words can not express the
gratitude I owe them.
College Park, April 2007
Zhenyu Zhang
iv
TABLE OF CONTENTS
List of Tables viii
List of Figures ix
1 Introduction 1
1.1 General overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Current state of the research . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2 Plasmon Resonances As an Eigenvalue Problem 14
2.1 Theory of (electrostatic) plasmon resonances . . . . . . . . . . . . . . 15
2.1.1 General considerations of (electrostatic) plasmon resonances . 15
2.1.2 Perturbation formulations . . . . . . . . . . . . . . . . . . . . 17
2.1.3 Zero order solution . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.4 First order correction . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.5 Second order correction . . . . . . . . . . . . . . . . . . . . . . 27
2.2 Comparison with the Mie theory for sphere . . . . . . . . . . . . . . . 29
2.2.1 Zero order solution . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.2 Second order solution . . . . . . . . . . . . . . . . . . . . . . . 32
2.3 Numerical analysis of plasmon resonances in nanoparticles . . . . . . 34
2.3.1 Numerical technique for solution of integral equation . . . . . 34
2.3.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.4 Plasmon resonances in shell structures . . . . . . . . . . . . . . . . . 47
2.4.1 A generalized eigenvalue problem . . . . . . . . . . . . . . . . 47
2.4.2 Analytical results for spherical and ellipsoidal shells . . . . . . 52
v
3 Fast Multipole Method for the Solution of Eigenvalue Problem 57
3.1 Why fast multipole method is needed? . . . . . . . . . . . . . . . . . 58
3.2 Factorization and translations of the kernel function . . . . . . . . . . 60
3.2.1 Factorization of the kernel function . . . . . . . . . . . . . . . 60
3.2.2 Translation in 3D . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.2.3 Difierentiation of the 1r kernel in 3D . . . . . . . . . . . . . . . 65
3.3 FMM algorithm and implementation . . . . . . . . . . . . . . . . . . 68
3.3.1 FMM algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.3.2 Date structure for FMM implementation . . . . . . . . . . . . 72
3.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.5 Multipole expansion method for spherical nanoparticles . . . . . . . . 82
3.5.1 The dual formulation for multiple spheres . . . . . . . . . . . 82
3.5.2 Multipole expansion method . . . . . . . . . . . . . . . . . . . 84
3.5.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . 92
4 Extinction Cross Sections of Nanoparticles 94
4.1 What is the extinction cross section ? . . . . . . . . . . . . . . . . . . 95
4.1.1 Concept of extinction cross section . . . . . . . . . . . . . . . 95
4.1.2 The optical theorem . . . . . . . . . . . . . . . . . . . . . . . 98
4.1.3 Measurement of ECS . . . . . . . . . . . . . . . . . . . . . . . 99
4.2 Computation of ECS . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.2.1 Near fleld computation . . . . . . . . . . . . . . . . . . . . . . 100
4.2.2 Far fleld computation . . . . . . . . . . . . . . . . . . . . . . . 106
4.3 Numerical techniques to compute ECS . . . . . . . . . . . . . . . . . 108
4.3.1 Discretization of the inhomogeneous integral equation . . . . . 108
vi
4.3.2 Far fleld and dipole moment . . . . . . . . . . . . . . . . . . . 109
4.3.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . 111
4.4 Numerical analysis of plasmon waveguides of light . . . . . . . . . . . 117
4.4.1 The numerical technique . . . . . . . . . . . . . . . . . . . . . 117
4.4.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . 120
5 Conclusion 126
Bibliography 130
vii
LIST OF TABLES
2.1 Eigenvalues for a single nanosphere. . . . . . . . . . . . . . . . . . . . 38
2.2 Eigenvalues for a single nanoellipsoid with aspect ratio 1:1:1.55. . . . 39
2.3 Comparison with experimental results [72] for gold nanoring. . . . . . 41
2.4 Comparison with experimental results [42, 73] for gold ellipsoidal
cylinder and triangular prism. . . . . . . . . . . . . . . . . . . . . . . 41
2.5 Comparison of the resonance wavelengths for Au spherical nanoshells
with silicon cores. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.6 Comparison of the resonance values of ?+ for spherical nanoshells,
?1 = 5?0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.7 Comparison of the resonance values of ?+ for ellipsoidal nanoshells,
?1 = 15:2?0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.1 Comparison of the resonance wavelengths for three silver spherical
nanoparticle system. . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.2 Computation times in seconds. . . . . . . . . . . . . . . . . . . . . . . 82
3.3 Computed eigenvalue for two spheres with truncation number p = 10. 91
3.4 Computed eigenvalue for two spheres with truncation number p = 10
using submatrix. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
3.5 ComparisonofnumericalresultsofBEMandMEMfortwo-nanosphere
and three-nanosphere conflgurations. . . . . . . . . . . . . . . . . . . 93
3.6 ComparisonofcomputationcostsofBEMandMEMfortwo-nanosphere
and three-nanosphere conflgurations. . . . . . . . . . . . . . . . . . . 93
4.1 Resonance frequencies for gold 7-sphere-chain waveguide. . . . . . . . 121
viii
LIST OF FIGURES
1.1 Darkfleld image of plasmon resonant silver nanoparticles. . . . . . . . 2
1.2 Schematic of plasmonic waveguide of light. . . . . . . . . . . . . . . . 3
1.3 Local plasmon photonic transistor proposed in [21]. . . . . . . . . . . 4
2.1 Measured relative permittivity [68] of silver as a function of wave-
length(frequency). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2 Measured relative permittivity [68] of gold as a function of wave-
length(frequency). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3 The dielectric nanoparticle bounded by surface S. . . . . . . . . . . . 17
2.4 Nanoparticles on substrate. The \dash particle" is the mirror image
of the actual particle on the substrate (see formula (2.25)). . . . . . . 22
2.5 Surface charge distribution (eigenfunction) for the flrst 10 eigenmodes
of single nanosphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.6 Surface charge distribution (eigenfunction) for the flrst 10 eigenmodes
of single nano-ellipsoidal. . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.7 Two nanospheres on substrate. . . . . . . . . . . . . . . . . . . . . . 42
2.8 Resonance wavelength for two nanospheres on substrate. . . . . . . . 42
2.9 Two ellipsoidal cylinders on substrate. . . . . . . . . . . . . . . . . . 43
2.10 Resonance wavelength for two ellipsoidal cylinders on substrate. . . . 43
2.11 The plasmon mode light intensity of flrst two resonance modes for
triangular nanoprism, colorbars indicate the ratio of optical near fleld
intensity normalized to incident flelds. . . . . . . . . . . . . . . . . . 44
2.12 Comparison of the plasmon resonance wavelength of (A) nanocube
ensemble in water and (B) single nanocube on a glass substrate. The
down arrows indicate the experimental results (432 nm and 452 nm
in (A), 395 nm and 457 nm for (B)) while the up arrows show the
computational results (420 nm and 478 nm in (A), 410 nm and 443
nm in (B)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.13 Examples of surface mesh used in computations for difierent geometry
of three dimensional nanoparticles. . . . . . . . . . . . . . . . . . . . 46
ix
2.14 (a) Schematic of a nanoshell with dielectric core; (b) the equivalent
problem for plasmon resonance in nanoshells. . . . . . . . . . . . . . 48
2.15 Plasmon resonant wavelengths for Au spherical nanoshells with dif-
ferent shell thickness. The core of these shells is silicon. . . . . . . . . 51
3.1 Space partitioning. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.2 Data structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3 Multipole expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.4 Local expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.5 Illustrationoflocal-to-local(RjR), multipole-to-local(SjR), andmultipole-
to-multipole (SjS) translations from expansion centerr1to expansion
center r2 (translation vector t = r2?r1). The star shows location
of a source. Expansion about center r1 is valid inside the lighter and
darker gray area, while expansion about center r2 is valid only in the
darker gray area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.6 Step 2 of FMM algorithm: Multipole to multipole translation. . . . . 69
3.7 Step 3 of FMM algorithm: Local to Local translation. . . . . . . . . . 70
3.8 Step 3 of FMM algorithm: Multipole to local translation. . . . . . . . 71
3.9 Step 3 of FMM algorithm: Multipole to local translation. . . . . . . . 71
3.10 A  ow chart of the standard FMM. . . . . . . . . . . . . . . . . . . . 72
3.11 Hierarchical cells at three levels in three dimensions. . . . . . . . . . . 73
3.12 The CPU time in seconds required for computation of problems of size
N using direct method and FMM. FMM computations using trunca-
tion number p = 5. Computation is performed on Dell workstation
with Pentium 3.0 GHz processor, 4.0 GB RAM. . . . . . . . . . . . . 78
3.13 Absolute error vs truncation numbers p2 for the case N = 14496.
Computation is performed on Dell workstation with Pentium 3.0 GHz
processor, 4.0 GB RAM. . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.14 CPU times vs group parameters s (number of points in the smallest
box) for N = 14496. Computation is performed on Dell workstation
with Pentium 3.0 GHz processor, 4.0 GB RAM. . . . . . . . . . . . . 80
3.15 Three-sphere system for the nano-lens proposed in [83]. . . . . . . . . 81
x
3.16 The schematic of a 7-sphere cluster. . . . . . . . . . . . . . . . . . . . 81
3.17 Multiple spherical nanoparticles located in proximity. . . . . . . . . . 82
3.18 Points x and y are located at same sphere. . . . . . . . . . . . . . . . 86
3.19 Points x and y are located at difierent spheres. . . . . . . . . . . . . 87
3.20 Truncation number as a function of the gap between two spheres. . . 90
4.1 Scattering diagram of a nanoparticle. . . . . . . . . . . . . . . . . . . 95
4.2 Schematic illustration for the measurement of extinction cross sections.100
4.3 Computational results compared with the Mie theory predictions for
the extinction cross section of a single Au nanoparticle (diameter =
10 nm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.4 Computational results compared with the Mie theory predictions for
the extinction cross section of a single Ag nanosphere (diameter = 20
nm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.5 Computational results of the extinction cross section of a single Ag
nanoellipsoid, with difierent aspect ratios: 1:1:0.8, 1:1:1. 1:1:1.2.
1:1:1.4. 1:1:1.6 for curve 1, 2, 3, 4, 5, respectively. . . . . . . . . . . . 114
4.6 Computed extinction cross sections of Au nanorings placed on sub-
strate, the dimensions of these rings are provided in Table 2.3. . . . . 115
4.7 Computational extinction cross section of a single InSb nanosphere
(diameter = 20 nm) placed on a glass substrate ? = 2:25?0. . . . . . . 116
4.8 Schematic of plasmonic waveguide of light. . . . . . . . . . . . . . . . 118
4.9 Computed extinction cross sections of gold 7-sphere-chain waveguide
for the transverse and longitudinal polarization modes (diameter =
50nm, center-to-center distance = 25nm). . . . . . . . . . . . . . . . . 122
4.10 Schematic of gold 11-sphere-chain waveguide of light . . . . . . . . . 123
4.11 Computed extinction cross sections of gold 11-sphere-chain waveguide
for the longitudinal polarization mode (diameter = 50nm, center-to-
center distance = 25nm). . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.12 Schematic of gold L-shape 11-sphere waveguide of light . . . . . . . . 124
xi
4.13 ComputedextinctioncrosssectionsofgoldL-shape11-spherewaveguide
for the longitudinal polarization mode (diameter = 50nm, center-to-
center distance = 25nm). . . . . . . . . . . . . . . . . . . . . . . . . . 124
4.14 Schematic of gold T-shape 11-sphere waveguide of light . . . . . . . . 125
4.15 ComputedextinctioncrosssectionsofgoldT-shape11-spherewaveguide
for the longitudinal polarization mode (diameter = 50nm, center-to-
center distance = 25nm). . . . . . . . . . . . . . . . . . . . . . . . . . 125
xii
Chapter 1
Introduction
1.1 General overview
Due to their reduction of size and dimensionality, nanoparticles exhibit opti-
cal properties that difier in fundamental and valuable ways from the properties of
individual atoms, molecules and bulk matter [1]. For example, the bright colors in
stained glass windows or colloidal solutions are caused by small noble metallic par-
ticles [2]. This optical efiect is due to the strong interaction between incident light
and conduction electrons conflned to the small volume of the metallic nanoparticle.
Under the in uence of an oscillating electric fleld, the negatively charged conduction
electrons perform a collective motion with respect to the positive-ion background,
creating an efiective charge at the surface that results in a restoring force. The
electron oscillations are therefore called surface plasmon resonances. The resonance
modes are localized near the surface of the nanoparticle and exponentially decay
along the depth from the surface. Surface plasmon resonances in metallic nanopar-
ticles lead to strong light scattering and absorption and remarkable enhancement of
local electromagnetic flelds [3].
The optical properties of small metallic particles have fascinated many re-
searchers in the past. For instance, Faraday flrst recognized that the red color in
stained glasses is due to small-sized Au particles in 1857, and this peculiar phenom-
1
Figure 1.1: Darkfleld image of plasmon resonant silver nanoparticles.
enon was theoretically explained by the groundbreaking work of Mie in 1908 [4].
Nevertheless, the research on surface plasmon resonances has increased in volume
considerably in the past decade. This recent revival of interest is caused by the
advances and new capabilities in fabrication methods to produce metallic nanos-
tructures of well-deflned sizes and shapes [5, 6, 7], and subsequently the possibilities
to utilize plasmon resonant nanoparticles in scientiflc and technological applications
ranging from optics and sensors to biology and medicine.
A key feature of plasmonic nanostructures is that light can be conflned and
manipulated on a scale below the difiraction limit governed by conventional optics
[8]. This opens opportunities for new optical devices and technology. For example,
plasmon-resonant metal nanoparticles can be used to construct miniature optical
waveguides. The plasmon waveguides consist of arrays (chains) of metallic nanopar-
ticles (see Figure 1.2) with their plasmon resonance frequencies in the region of
2


	




	

Figure 1.2: Schematic of plasmonic waveguide of light.
optical waveguiding. At speciflc frequencies, incident optical radiation can excite
plasmon resonances in the flrst nanoparticle, which then through near-fleld coupling
can induce plasmon resonances in all nanoparticles that form the chain [9]. Atwa-
ter et al have demonstrated that optical signals can propagate over a micrometer
along the waveguide [10, 11, 12]. Surface plasmon resonances can also be exploited
in the fleld of optical lithography for the fabrication of nanoscale features beyond
the difiraction limit [13, 14, 15]. By using a photoresist that is sensitive at the
surface plasmon resonant frequency, the exposure of a thin layer of photoresist di-
rectly below a contact mask can create an aerial image on nanometer length scales.
This surface plasmon resonant nanolithography technique is not difiraction limited,
and can produce subwavelength features using broad beam illumination with visible
light.
Another salient feature associated with plasmon resonances is the strong local
fleld enhancement. For instance, plasmon-resonant metallic nanoparticles have been
used in surface enhanced Raman spectroscopy [16, 17, 18]. Moreover, active pho-
tonic devices have been proposed which utilize plasmon resonances. Mayergoyz and
3
Fredkin [19, 20] suggested that the optical controllability of semiconductor nanopar-
ticles can be utilized for the development of nanoscale light switches and all-optical
nanotransistors. On the other hand, Tominnaga and co-workers [21] showed that by
focusing two laser beams (405 and 635 nm) in one small spot on a high-speed rotat-
ing optical disk, a large signal enhancement is observed (Figure 1.3). It was found
that a plasmon interaction generated between a silver light-scattering nanoparticle
and small marks recorded in the optical disk with a super-resolution near-fleld struc-
ture produced the large signal ampliflcation in the spot (< 1 ?m). A modulated
signal of the blue laser was enhanced by 60 times by controlling the red laser power
from 1.5 to 3.5 mW. It has been shown that the system has the potential to realize
all-thin-fllms photonic transistors by using local plasmon ampliflcation.
Figure 1.3: Local plasmon photonic transistor proposed in [21].
Surface plasmon resonances are very sensitive to the local dielectric environ-
4
ment. In other words, a small change of the refractive index of surrounding medium
will result in appreciable shift of plasmon resonant frequencies as well as the extinc-
tion coe?cients. A variety of nanoscale biosensors have been demonstrated based on
this principle [22, 23, 24]. For example, Storhofi and Mirkin linked a single-stranded
DNA to a gold nanoparticle of 15 nm in diameter. When this DNA hybridizes with
its complementary DNA in the test sample, the duplex formation leads to aggre-
gation of nanoparticles, shifting the plasmon resonance from burgundy-red color
(single nanoparticle) to blue black (aggregation of nanoparticles) [25, 26]. In a
similar manner, Halas and co-workers used gold nanoshells to demonstrate a rapid
immunoassay, capable of detecting analyte within a complex biological media such
as blood [27, 28].
Another area which attracts considerable interest lately is the design and engi-
neering of metamaterials. Metamaterials are artiflcial composite materials that have
negative permeability and permittivity simultaneously, also known as left-handed
or double negative materials [29]. Due to their unique electromagnetic properties,
metamaterials are promising for many optical and microwave applications, such as
perfect lenses, transmission lines, and antennas. It is well known that many metals
(e.g. silver and gold) have negative ? at visible wavelengths. Relying on the plas-
mon resonances of individual metallic nanoparticles and arranging the geometry of
composite inclusions to resemble magnetic nanoloops, negative ? and left-handed
metamaterials with a negative index of refraction in the optical frequency range can
be realized [30, 31].
In addition to the aforementioned applications, the surface plasmon resonance
5
also has potential in next generation optical storage [32, 33, 34], scanning near-fleld
optical microscopy [35, 36] and energy solar cells [37, 38]. Since the understand-
ing of the optical properties of metallic nanoparticles holds both fundamental and
practical signiflcance, it is important to have the capability to model and predict
the plasmon resonance behavior of nanoparticles with complex shape. The goal of
this dissertation is to present a new and e?cient method for the analysis of plasmon
resonances in nanoparticles.
1.2 Current state of the research
In the classical picture, the optical property of nanoparticles is treated as a
scattering problem, which can be described by Maxwell?s equations. In this frame-
work, thecompletedescriptionofthematerialpropertiesisincludedinthedispersion
relation, which gives the complex permittivity ?(!) as a function of the frequency
or wavelength. At speciflc negative permittivity values, plasmon resonances will be
excited in these small particles. These speciflc permittivity values strongly depend
on the particle shape, size and surrounding medium, since the boundary conditions
imposed by Maxwell?s equations determine whether such a particle resonance can
build up. One may wonder whether the classical description of the material based
solely on Maxwell?s equations and bulk dispersion relations is appropriate for the
small particles. In fact, it has been experimentally demonstrated that this macro-
scopic approach is adequate for particle dimensions as small as a few nanometers
[3, 39]. Furthermore, this question has been discussed in a recent paper by Levi
6
et al [40]. According to their quantum mechanical computations, it can be con-
cluded that there is no appreciable discrepancy between quantum model and classic
electrodynamics for particle size above 8 nanometers.
The optical properties of small spherical metallic particles accounting for the
surface plasmon resonance were flrst explained theoretically by the groundbreaking
work of Mie in 1908 [4]. Mie solved Maxwell?s equations for an electromagnetic light
wave interacting with a small sphere that has the same macroscopic, frequency-
dependent material dielectric constant as the bulk metal. Mie?s electrodynamics
computation gave a series of multipole oscillations for the extinction and scattering
cross sections of the particles as a function of the particle radius. For small but
flnite-sized particles, Mie theory gives the resonance values of dielectric permittivity
of the spherical particle as [39]:
?+ = ?
 
2+ 125 !2?0?0a2
?
?0; (1.1)
where ! is the angular frequency of the optical radiation, a is the radius of the
sphere, and ?0 and ?0 have their usual meaning. Mie?s results are still widely used
today since it is the only available analytical solution and most of the nanoparticles
produced by the colloidal or wet chemistry method are more or less spherical.
As soon as one considers particles other than spheres or even spheres within
an asymmetric dielectric environment, it is usually not possible to obtain analytical
solutions to Maxwell?s equations. Due to this fact, there has been a great deal of
efiort put into developing numerical methods. Indeed, the literature on numerical
methods for classical electrodynamics problem is enormous. However, there are some
7
important challenges and simpliflcations associated with metallic nanoparticles that
limit the applicability of methods used to study problems with much longer wave-
length and length scales. In recent years, there are several numerical techniques that
have been used to describe optical properties of non-spherical metallic nanoparticles.
They can be brie y summarized as follows.
(1) Discrete dipole approximation (DDA) method. Schatz and co-workers
[41, 42] have applied the DDA method to non-spherical nanoparticles which can be
traced to Purcell [43] and Draine [44]. In DDA, a particle with arbitrary shape is
treated as a three-dimensional assembly of dipoles on a cubic grid. Each dipole cell
is assigned a complex polarizability fii which can be computed from the complex
dispersion relation and the number of dipoles in a unit volume. This polarizability
causes an oscillating dipole moment or a polarization Pi at each cell, depending on
the total electric fleld at the respective position and given by
Pi = fii ?Ei: (1.2)
The total fleld Ei is the sum of the incident fleld Einc;i (usually a plane wave) and
a contribution from all other dipoles Edip;i:
Ei = Einc;i +Edip;i = E0eik?ri +X
i6=j
Aij ?Pj: (1.3)
The matrix Aij includes the interaction of all dipoles depending on the distance
of the dipoles; thus, a full dense matrix results. This equation can be solved by
iteration. Once the polarizations Pi are found, the extinction and absorption cross
sections can be computed by using standard formulas. The DDA technique has
been applied to difierent shape of nanoparticles and was able to reproduce the
8
spectral shape of the plasmon resonances. It has two weaknesses, however. The
total volume of material that can be described is limited by available computer
resources to dimensions of a few hundred nanometers and the electric flelds close to
particle surface are inaccurate [41].
(2) Finite-difierence time-domain (FDTD) methods. The popular FDTD tech-
nique has also been used for the modeling of plasmon resonances [45, 46, 47]. In this
technique, the Maxwell?s equations in difierential form are solved directly. The space
and time are both discretized and a difierence approximation is applied to evaluate
the space and time derivatives of the fleld. The equations are solved in a leap-frog
manner; that is, the electric fleld is solved at a given instant in time, then the mag-
netic fleld is solved at the next instant in time, and the process is repeated over and
over again. Although the FDTD technique is rather versatile and easy to implement,
it has several drawbacks as far as as the analysis of plasmon resonance modes in
nanoparticles is concerned. First, it is a \probing" technique when a nanoparticle
is \numerically probed" many times by incident radiation of difierent frequencies
and polarizations. In other words, FDTD is a \trial-and-error" approach to the
analysis of plasmon resonance modes in nanoparticles, and this is detrimental to the
efiectiveness of FDTD for plasmon mode computations. Second, FDTD requires dis-
cretization of three-dimensional space, since only a flnite region of three-dimensional
space can be discretized and used in computations, FDTD requires the introduction
of artiflcial external boundaries and special absorbing boundary conditions on these
boundaries to minimize distortions and errors caused by the introduction of artiflcial
boundaries. Finally, plasmon resonances occur in dispersive dielectric media with
9
non-instantaneous in time (convolution-type) constitutive relations between electric
displacement and electric fleld. These non-instantaneous constitutive relations lead
to flnite difierence schemes with the electric fleld coupling at all previous time-steps.
This past history coupling of electric fleld diminishes the efiectiveness of FDTD for
plasmon resonance computations.
(3) T-matrix method [48, 49, 50]. In T-matrix method, the incident and
scattered electric flelds are expended in series of suitable vector spherical wave func-
tions, and the relation between the columns of the respective expansion coe?cients
is established by means of a transition matrix (or T-matrix). The elements of the
T-matrix are independent of the incident and scattered electric flelds and only de-
pend on the shape, size, refractive index of the particle as well as its orientation
with respect to the reference frame. These elements are obtained by numerical
integration. While the T-matrix method is very accurate and powerful, it sufiers
from two major disadvantages. First, for a particle with arbitrary shape, the el-
ements of T-matrix have to be computed through a surface integral. As this is
computationally expensive, most implementations of this method are restricted to
axisymmetric scatterers. In this case, the surface integrals can be reduced to line
integrals. Second, the numerical stability of this method is compromised for highly
aspherical particles, for which the convergent size of the T-matrix should be large,
and T-matrix computations may converge very slowly or even diverge.
(4) Multiple multipole (MMP) method [51, 52, 53]. In MMP method, the flelds
in the individual domains of the structure under investigation are described by series
expansion of known analytical solutions of Maxwell equations (multipole functions,
10
plane waves, waveguide modes, etc.). Each solution has a free parameter that can
be determined from the boundary conditions to be satisfled in a suitably selected
set of points on the interfaces between the domains. Usually more conditions than
free parameters are used, and the resulting over-determined system of equations is
solved in the least-squares sense. With this procedure a smooth error distribution is
obtained on the boundaries, and one can evaluate the errors in the individual bound-
ary points to estimate the quality of the solution. Note that Maxwell equations are
exactly fulfllled inside the domains but are only approximated on the boundaries.
In this case, the computational efiort varies with the number of matching points,
expansion functions, and basis functions used in the expansions. The computational
advantage of the MMP is that only the boundaries need to be discretized and not
the domains themselves.
It can be seen that by using the above mentioned methods, plasmon resonances
in metallic nanoparticles are found numerically by using a \trial-and-error" method,
i.e., by probing metallic nanoparticles of complex shapes with radiation of various
frequencies and polarizations. A mode-based approach and the direct calculations
of plasmon resonance modes are clearly preferable. In this dissertation, an integral
equation based numerical technique is presented to fully characterize the plasmon
resonances in nanoparticles [9, 20, 58, 59, 60, 61, 62, 63]. The basic idea of the tech-
nique is to frame the plasmon resonance as an integral equation based eigenvalue
problem, which was proposed by Mayergoyz and Fredkin [19] and can be traced
back to publications of Mayergoyz [54, 55]. The integral equation technique has
been also used for the analysis of plasmon resonances in a limited form by Ouyang
11
and Isaakson [56, 57]. By using this technique, the plasmon resonance modes and
resonant frequencies of three dimensional nanoparticles can be computed directly
at the quasi-electrostatic limit by solving an eigenvalue problem for speciflc homo-
geneous surface integral equations. The calculation of high order corrections which
counts for the flnite wavelength nature of radiation has been developed through
perturbation technique. The integral equation technique has been also extended
to the computation of extinction cross-sections of nanoparticles, where the problem
is reduced to the solution of an inhomogeneous integral equation. Furthermore,
in this dissertation work, the fast multipole method (FMM) [64, 65] is applied for
the large scale eigenvalue problem to reduce the computational cost from O(N2) to
O(NlogN) which is a great advance when N is large. For these reasons, the integral
equation based numerical technique is very e?cient and can be instrumental for the
design of plasmon resonant nanoparticles and to tailor their optical properties for
various applications.
1.3 Outline
This dissertation is organized as follows. Chapter 2 presents the basic idea
of the electrostatic (plasmon) resonance (ESR) model. It is shown that the plas-
mon resonant frequency of a nanoparticle can be directly computed through the
solution of an integral equation based eigenvalue problem. The technique for the
discretization of the integral equation is discussed and numerous numerical results
are presented and compared with theoretical results as well as available experimental
12
data for various three-dimensional nanoparticles. In the last part of this chapter, the
ESR model is extended for the analysis of plasmon resonances in metallic nanoshells,
which is very promising for biological sensing applications.
The numerical aspect of computations of the eigenvalue problem is consid-
ered in Chapter 3. The fast multiple method (FMM) is introduced for the large
scale eigenvalue problem. The applicability of FMM and its implementation in
the iterative algorithm are presented in detail. Numerical results which illustrate
the e?ciency of FMM and a comparison with traditional methods for the numerical
analysis of plasmon resonances are presented. In addition, a semi-analytical method
based on multipole expansions is developed for the analysis of plasmon resonances
in spherical nanoparticles.
In Chapter 4, the computation of extinction cross sections (ECS) of nanopar-
ticles is demonstrated. In this chapter, the plasmon resonance in nanoparticles is
treated as a scattering problem. It is shown that the scattered electromagnetic
flelds can be obtained via the solution of an inhomogeneous integral equation and
then the ECS is computed by invoking the optical theorem. Simulation results are
presented and compared with theoretical results and available experimental data
for three-dimensional nanoparticles. In the last section of Chapter 4, the plasmonic
waveguide of light is discussed and numerically investigated.
Finally, conclusions and further work are presented in Chapter 5.
13
Chapter 2
Plasmon Resonances As an Eigenvalue Problem
In this chapter, the theory and numerical analysis of (electrostatic) plasmon reso-
nance modes in three-dimensional nanoparticles is presented. It is shown that the
problem of determining plasmon resonance modes in nanoparticles can be treated as
an integral-equation-based classical eigenvalue problem, and the plasmon resonant
frequencies can be directly found through the solution of this eigenvalue problem
(sections 2.1). It is also demonstrated that radiation corrections due to flnite sizes
of nanoparticles can be computed through perturbation technique. The comparison
between this method and the Mie theory for spherical nanoparticle is demonstrated
in section 2.2. The numerical technique to solve the integral equation as well as
the numerical analysis of plasmon resonances in three-dimensional nanoparticles are
presented in section 2.3. Finally, the integral equation technique is extended to the
generalized eigenvalue problem for the analysis of nanoshell structures in section
2.4.
14
2.1 Theory of (electrostatic) plasmon resonances
2.1.1 General considerations of (electrostatic) plasmon resonances
It is known that plasmon resonances in nanoparticles occur at speciflc frequen-
cies for which the particle permittivity is negative and the free-space wavelength of
the radiation is large in comparison with particle dimensions. The latter condition
suggests that these resonances are electrostatic in nature. Indeed, since particle
dimensions are small compared to the free-space wavelength of the radiation (when
resonances occur), time-harmonic electromagnetic flelds within the nanoparticles
and around them vary almost with the same phase. As a result, at any instant
of time these flelds look like electrostatic flelds. Therefore, in the study of plas-
mon resonances in nanoparticles, we shall follow the traditional approach where all
losses are flrst neglected and resonance frequencies are found for lossless systems as
frequencies for which source-free electromagnetic flelds may exist. This approach
leads to the consideration of resonances in the electrostatic limit where all radiation
losses are flrst neglected. When the dielectric permittivity of nanoparticles is nega-
tive, the uniqueness theorem of electrostatics is broken. For this reason, source-free
electrostatic flelds may appear for certain negative values of dielectric permittivi-
ties, which is the manifestation of resonances. The frequencies corresponding to the
above negative values of permittivity are the resonance frequencies.
It is clear from the previous discussion that electrostatic resonances may occur
only in particles whose media exhibit dispersion and the real part of the permittivity
assumes negative value for some range of frequencies. For metals, this frequency
15
0 200 400 600 800 100012001400160018002000?200
?150
?100
?50
0
50
Wavelength (nm)
Permittivity of Silver
Real partImaginary part
Figure 2.1: Measured relative permittivity [68] of silver as a function of wave-
length(frequency).
0 200 400 600 800 100012001400160018002000?200
?150
?100
?50
0
50
Wavelength (nm)
Permittivity of Gold
Real partImaginary part
Figure 2.2: Measured relative permittivity [68] of gold as a function of wave-
length(frequency).
16
range is below the plasma frequencies !p, but at su?ciently high frequencies for
collisions to be unimportant. The dispersion relation of many metals (but not all)
can be well described by the Drude model:
?(!) = ?0
"
1? !
2
p
!(! +i )
#
; (2.1)
!2p = nee
2
?0me: (2.2)
where  is the damping factor or electron relaxation rate. The most commonly used
data from literature for ?(!) are the tables in Paliks Handbook of Optical Constants
[69] and the tables by Johnson and Christy [68]. As an example, the dispersion
relations for silver and gold [68] are shown in Figures 2.1 and 2.2, respectively.
2.1.2 Perturbation formulations
()??? ++= 0??=
?+V
?VS
n

?
Figure 2.3: The dielectric nanoparticle bounded by surface S.
Consider a nanoparticle of arbitrary shape with uniform dielectric permittivity
?+(!) (Fig. 2.3). We are interested in negative values of ?+(!) for which a source-
free electromagnetic fleld may exist. To flnd such permittivities, we shall start with
17
Maxwell equations in term of the vectors (where e?i!t is assumed for E, H, e, h):
e = ?120E; h = ?120H (2.3)
and spatial coordinates scaled by the general diameter d of the nanoparticle. This
leads to the following boundary value problem:
r?e+ = ?iflh+; r?h+ = i?+?
0
fle+; (2.4)
r?e+ = 0; r?h+ = 0; (2.5)
r?e? = ?iflh?; r?h? = ifle?; (2.6)
r?e? = 0; r?h? = 0; (2.7)
n?
?
e+ ?e?
?
= 0; n?
?
h+ ?h?
?
= 0; (2.8)
n?
 ?
+
?0 e
+ ?e?
?
= 0; n?
?
h+ ?h?
?
= 0; (2.9)
where superscripts \+" and \-" are used for physical quantities inside (V +) and
outside (V ?) the nanoparticle, respectively, n is the outward unit normal to S, and
fl = !p?0?0d: (2.10)
Since the free-space wavelength ? is large in comparison with the dimension of
nanoparticle, fl can be treated as a small parameter and source-free solutions of
the boundary value problem (2.4){(2.9) and dielectric permittivities ?+(!) at which
they occur can be expanded in terms of fl:
e? = e?0 +fle?1 +fl2e?2 +:::; (2.11)
h? = h?0 +flh?1 +fl2h?2 +:::; (2.12)
18
?+ = ?(0)+ +fl?(1)+ +fl2?(2)+ +:::: (2.13)
By substituting formulas (2.11){(2.13) into equations (2.4){(2.7) as well as boundary
conditions (2.8){(2.9) and equating terms of equal powers of fl, we can obtain the
boundary value problems for e?k and h?k (k = 0;1;2;:::).
2.1.3 Zero order solution
For zero-order terms, these boundary value problems are written in terms of
E?0 = ??120 e?0 and H?0 = ??120 h?0 as follows:
r?E?0 = 0; r?E?0 = 0; (2.14)
n?
?
E+0 ?E?0
?
= 0; n?
0
@?
(0)
+
?0 E
+
0 ?E
?
0
1
A = 0; (2.15)
and
r?H?0 = 0; r?H?0 = 0; (2.16)
n?
?
H+0 ?H?0
?
= 0; n?
?
H+0 ?H?0
?
= 0: (2.17)
From equations (2.16) and (2.17), it is clear that:
H?0 = 0: (2.18)
To solve E0 from (2.14)-(2.15), an electric potential ? can be introduced for
this source-free electric fleld and the potential can be represented as the electric
potential of single layer of electric charges  distributed over the boundary S of the
particle (Fig. 2.3):
?(Q) = 14??
0
I
S
 (M)
rMQ dSM; (2.19)
19
where Q and M are the observation point and integration point on boundary S,
respectively. In other words, a single layer of electric charges on S may create the
same electric fleld in the free space as the source-free electric fleld that may exist
in the presence of the dielectric nanoparticle with negative permittivity. It is clear
that the electric fleld of surface charges  is curl and divergence free in V + and
V ? and satisfles the flrst boundary condition in (2.15). Next, we recall that the
normal components of electric fleld induced by surface electric charges are given by
the formulas[70, 71]:
n(Q)?E?0 (Q) = currency1 (Q)2?
0
+ 14??
0
I
S
 (M) rMQ ?nQr3
MQ
dSM: (2.20)
By substituting formulas (2.20) into the second boundary condition (2.15), after
simple transformations we arrive at the following homogeneous boundary integral
equation:
 (Q) = ?2?
I
S
 (M) rMQ ?nQr3
MQ
dSM; (2.21)
where the eigenvalue ? (not wavelength) is deflned by:
? = ?
(0)
+ (!)??0
?(0)+ (!)+?0: (2.22)
Thus, source-free electric flelds may exist only for such values of permittivity ?(0)+ that
the integral equation (2.21) has nonzero solutions. In other words, in order to flnd
the resonance values of ?(0)+ (and the corresponding resonance frequencies) as well as
resonance modes, the eigenvalues and eigenfunctions of the integral equation (2.21)
must be computed. Integral equation (2.21) and its spectrum for the calculations
of electrostatic and scattering problem with negative ? were extensively studied in
20
publications of Mayergoyz [54, 55] and later introduced for the analysis of plasmon
resonances by Ouyang and Isaacson [56, 57]. Then it was greatly expanded in
publications [19, 9, 20, 58, 59, 60, 61, 62, 63]
For particles of complex shapes the resonance frequencies and resonance modes
can be found through numerical solution of integral equation (2.21). If the boundary
S of the particle is not smooth, then  (M) may have singularities at the corners and
the edges of S that may negatively afiect the accuracy of numerical computations.
In this situation, the dual formulation can be employed, where dipole (double layer
of electric charges) density ? (M) is distributed over S in such a way that it creates
the same electric displacement fleld (D0) in free space as the source-free electric
displacement fleld that may exist in the presence of the dielectric particle with neg-
ative permittivity. By using the known properties of double-layer potential[70, 71],
it can be shown that the two displacement flelds mentioned above will be identical
if the dipole density ? (M) satisfles the following homogeneous boundary integral
equation:
? (Q) = ?2?
I
S
? (M) rQM ?nMr3
QM
dSM; (2.23)
where ? is given by formula (2.22). It is apparent that the boundary integral equa-
tion (2.23) is adjoint to the integral equation (2.21). For this reason, it has the
same spectrum as one can expect on the physical grounds. The dipole density ?(M)
is proportional to the discontinuity of double-layer potential across S and, conse-
quently, it is flnite even for non-smooth boundaries S. This is the advantage of
integral equation (2.23) for numerical computations with non-smooth boundaries
21
S?
substrate
air
M?
M
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a0 a0 a0 a0 a0 a0 a0 a0 a0 a0
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
a1 a1 a1 a1 a1 a1 a1 a1 a1 a1
Figure 2.4: Nanoparticles on substrate. The \dash particle" is the mirror image of
the actual particle on the substrate (see formula (2.25)).
S. In applications, nanoparticles are usually placed on dielectric substrates (see
Fig. 2.4). In this case, the integral equation (2.21) can be modifled as follows:
 (Q) = ?2?
I
S
 (M)n(Q)?rQ [G(Q;M)]dSM; (2.24)
where ? is given by formula (2.22), while G(M;Q) is the Green function deflned by
the formula:
G(Q;M) = 1r
MQ
? ???0?+?
0
1
rM0Q (2.25)
Here, ? is the permittivity of the substrate and M0 is the image of M with respect
to the substrate plane eS.
The presented discussion can be easily extended to the analysis of electrostatic
(plasmon) resonances of several particles located in proximity to one another. In
this case S in integral equations (2.21) and (2.23) must be construed as the union
of boundaries of all dielectric particles, while  (M) (and ? (M)) are deflned on this
union. It has to be noted, however, as the number of nanoparticles increases, the
computational cost for solving the integral equations (2.21) and (2.23) becomes very
expensive and this bottleneck problem will be addressed in Chapter 3.
22
The kernels in integral equations (2.21) and (2.23) have weak (integrable in
the usual sense) singularities. For this reason, the integral operators in the above
integral equations are compact. This implies that the plasmon spectrum is discrete
despite the inflnite region of fleld distribution. It has been proved [70, 71] that the
spectrum of integral equation (2.21) has the following properties: for any shapes
of S all eigenvalues are real, ? = 1 is an eigenvalue, and for all other eigenvalues
j?j > 1. It is apparent from (2.21) that the eigenvalue ? = 1 corresponds to the
case of ?(0)+ ! 1, and the respective eigenfunction  (M) can be construed as the
distribution of surface electric charges over the boundary S of a conductor V +. This
eigenvalue is irrelevant as far as the discussion of electrostatic (plasmon) resonances
is concerned. All other eigenvalues correspond to source-free (resonance) conflgura-
tions of electrostatic flelds and, according to (2.21), these conflgurations may exist
(as expected) only for negative values of ?(0)+ . After these negative resonance values
of ?(0)+ are found through the solution of integral equation (2.21), the appropriate
dispersion relation can be employed to flnd the corresponding resonance frequencies.
When losses are not neglected, actual permittivities are complex-valued functions
of frequency !. These permittivities may assume real resonance values ?(0)+ only for
complex resonance frequencies.
The kernels in integral equations (2.21) is not Hermitian (not self-adjoint),
because the kernel of this equation is not symmetric. For this reason, the eigen-
functions  i(M) and  k(M) corresponding to difierent eigenvalues ?i and ?k are not
orthogonal on S in the usual sense. Nevertheless, it has been shown [20] that elec-
trical flelds E0i and E0k corresponding to eigenfunctions  i(M) and  k(M) satisfy
23
the following strong orthogonality conditions:
Z
V?
E0i ?E0kdV = 0: (2.26)
Furthermore, the eigenfunctions  k(M) and ?i(M) of adjoint equations 2.21) and
(2.23) form two biorthogonal sets
I
S
 k(M)?i(M)dS = ?ki: (2.27)
This biorthogonal property is very important and can be used for the analysis of
time-dynamics of plasmon resonances in nanoparticles [61].
It is instructive to note that the \plasmon" eigenfunctions of integral equation
(2.21) have the property:
I
S
 (M)dSM = 0: (2.28)
Indeed, by integrating both sides of equation (2.21) with respect to Q and by using
the facts that:
I
S
rMQ ?nQ
r3MQ dSM = 2? (2.29)
and ? 6= 1 for plasmon resonances, we arrive at (2.28).
It is apparent that the mathematical structure of integral equation (2.21) is
invariant with respect to the scaling of S, i.e. the scaling of the dimensions of the
particle. This leads to the unique property of electrostatic (plasmon) resonances:
resonance frequencies depend on particle shape but they are scale invariant with
respect to particle dimensions, provided that they remain appreciably smaller than
the free-space wavelength.
24
2.1.4 First order correction
From formulas (2.4){(2.9), (2.11){(2.13) and (2.18), we derive the following
boundary value problems for the flrst-order corrections e?1 , ?(1)+ and h?1 , respectively:
r?e?1 = 0; r?e?1 = 0; (2.30)
n?
?
e+1 ?e?1
?
= 0; n?
0
@?
(0)
+
?0 e
+
1 ?e
?
1
1
A = ??
(1)
+
?0 n?e
+
0 ; (2.31)
and
r?h+1 = i?
(0)
+
?0 e
+
0 ; r?h
?
1 = ie
?
0 ; r?h
?
1 = 0; (2.32)
n?
?
h+1 ?h?1
?
= 0; n?
?
h+1 ?h?1
?
= 0; (2.33)
where as before e?0 = ?120E?0 .
To solve (2.30)-(2.31), the electric potential ?1 of a single layer of electric
charges  1 distributed over S can be introduced for the electric fleld e?1 . Then,
by using formulas (2.20) and the same reasoning as in the derivation of integral
equation (2.21), we arrive at the following integral equation for  1 (M):
 1 (Q)? ?2?
I
S
 1 (M) rMQ ?nQr3
MQ
dSM = ?(1)+ 2?0?(0)
+ +?0
nQ ?e+0 (Q); (2.34)
where ? is given by formula (2.22).
It is clear that ? is one of the eigenvalues of integral equation (2.21), because
only for such ? nonzero fleld e0 exists. Since ? in (2.34) is an eigenvalue, a solu-
tion to equation (2.34) exists only under the condition that the right-hand side of
equation (2.34) is orthogonal on S to a nonzero solution ? (Q) of the correspond-
ing homogeneous adjoint equation (2.23) with the same eigenvalue ? (this is the
25
so-called \normal solvability condition"). It is clear that nQ?e+0 (Q) is proportional
to @?+@n (Q). By using the well-known properties of double-layer potential [70, 71], it
can be shown that ? (Q) is proportional to ?+ (Q). Consequently,
I
S
? (Q)nQ ?e+0 dSQ = fi
I
S
?+ (Q) @?
+
@n (Q)dSQ = fi
I
V +
jr?+j2dV 6= 0: (2.35)
This means that the integral equation (2.34) is only solvable if
?(1)+ = 0: (2.36)
Thus, for any shape of nanoparticle the flrst order correction for resonant values
of dielectric permittivity is equal to zero. As a result of (2.36), integral equation
(2.34) is reduced to homogeneous integral equation identical to equation (2.21). This
implies that up to a scale  (M) and  1 (M) as well as e?0 and e?1 are identical. For
this reason, it can be assumed that:
e?1 = 0: (2.37)
Next, we proceed to the solution of boundary value problem (2.32){(2.33).
Terms i?
(0)
+
?0 e
+
0 and ie
?
0 in the flrst two equations (2.32) can be interpreted as current
sources and the solution of boundary value problem (2.32){(2.33) can be written in
the following integral form:
h1 (Q) = i?
(0)
+
4??0
Z
V +
e+0 (M)?rMQ
r3MQ dVM +
i
4?
Z
V?
e?0 (M)?rMQ
r3MQ dVM: (2.38)
The last expression can be appreciably simplifled and reduced to an integral over
the boundary S. Indeed, by using the fact that r?e?0 = 0 and by employing the
\curl-theorem"[87], after simple transformations we arrive at:
h1 (Q) = ?
i
 
?(0)+
?0 ?1
?
4?
I
S
nM ?e0 (M)
rMQ dSM: (2.39)
26
2.1.5 Second order correction
In this section, we proceed to the discussion of second order corrections for
?+ (!). From formulas (2.4){(2.9), (2.11){(2.13) and (2.36){(2.37) we derive the
following boundary value problems for the second-order corrections of e?2 , ?(2)+ and
h?2 , respectively:
r?e?2 = ?ih?1 ; r?e?2 = 0; (2.40)
n?
?
e+2 ?e?2
?
= 0; n?
0
@?
(0)
+
?0 e
+
2 ?e
?
2
1
A = ??
(2)
+
?0 n?e
+
0 ; (2.41)
and
r?h?2 = 0; r?h?2 = 0; (2.42)
n?
?
h+2 ?h?2
?
= 0; n?
?
h+2 ?h?2
?
= 0; (2.43)
It is apparent from (2.42){(2.43) that
h?2 = 0: (2.44)
To solve (2.40){(2.41), we shall flrst decompose the electric fleld e?2 into two distinct
components:
e?2 = ~e?2 + ~~e?2 ; (2.45)
which satisfy the following boundary value problems, respectively:
r?~e2 = ?ih1; r?~e2 = 0; (2.46)
n?
?
~e+2 ?~e?2
?
= 0; n?
?
~e+2 ?~e?2
?
= 0; (2.47)
and
r?~~e?2 = 0; r?~~e?2 = 0; (2.48)
27
n?
?~
~e+2 ?~~e?2
?
= 0;
n?
 
?(0)+
?0 ~~e
+
2 ?~~e
?
2
?
= ?n?
?
?(2)+
?0 e
+
0 +
 
?(0)+
?0 ?1
?
~e+2
?
: (2.49)
The solution of the boundary value problem (2.46){(2.47) is:
~e2 (P) = ? i4?
Z
R3
h1 (Q)?rQP
r3QP dVQ: (2.50)
By substituting the expression (2.39) into the last formula and by changing the
order of integration, we arrive at:
~e2 (P) = ?
 
?(0)+
?0 ?1
?
16?2
I
S
(nM ?e0 (M)) ?
?Z
R3
rQP
r3QP
1
rQM dVQ
!
dSM: (2.51)
The last expression can be simplifled to [20]:
~e2 (P) = ?
 
?(0)+
?0 ?1
?
8?
I
S
(nM ?e0 (M))?rMP
rMP dSM: (2.52)
Next, we proceed to the solution of the boundary value problem (2.48){(2.49). It
is apparent that the electric potential ?2 of single layer of electric charges  2 (M)
distributed over S can be introduced for the electric fleld ~~e?2 . Then, by using formula
(2.20) and the same line of reasoning as in the derivation of integral equation (2.21),
we obtain the following integral equation for  2 (M):
 2 (Q)? ?2?
I
S
 2 (M) rMQ ?nQr3
MQ
dSM = 2?
2
0
?0 ??(0)+ n?
2
4?
(2)
+
?0 e
+
0 +
0
@?
(0)
+
?0 ?1
1
A~e+2
3
5;
(2.53)
where as before ? is given by formula (2.22) and is one of the eigenvalues of inte-
gral equation (2.21). Since ? is an eigenvalue, a solution to equation (2.53) exists
(according to the Fredholm Theorem) only under the condition that the right-hand
side of the equation (2.53) is orthogonal on S to a nonzero solution ? (Q) of the
28
corresponding homogeneous adjoint equation (2.23). This normal solvability con-
dition leads to the following expression for the second order correction of resonant
permittivity ?(2)+ :
?(2)+ =
?0
 
?(0)+
?0 ?1
?H
S ? (Q)nQ ?~e
+
2 (Q)dSQ
H
S ? (Q)nQ ?e
+
0 (Q)dSQ
: (2.54)
Thus, the algorithm for computations of the second order correction ?(2)+ can be
stated as follows: flrst, integral equations (2.21) and (2.23) are solved and for each
eigenvalue the corresponding ?(0)+ and eigenfunction  (M) and ? (M) are found;
then, by using formula (2.52), ~e+2 is computed on S; flnally, by employing expres-
sion (2.54), the second-order corrections ?(2)+ to resonant permittivities ?(0)+ can be
calculated. According to (2.13), the resonant permittivities are given by the formula:
?+ = ?(0)+ +fl2?(2)+ ; (2.55)
where fl is deflned in (2.10).
2.2 Comparison with the Mie theory for sphere
The perturbation formulations have been derived in section 2.1 for the cal-
culation of the values of dielectric permittivity ?+ for speciflc plasmon modes. To
test its validity, let us flrst compare with the Mie theory, which provides analytical
solutions for spherical nanoparticles. Indeed, according to the Mie theory, the res-
onance values of dielectric permittivity (up to the second order corrections) for the
flrst three degenerate plasmon modes (spatially uniform modes) is [39]:
?+ = ?
 
2+ 35fl2
?
?0 = ?
 
2+ 125 !2?0?0a2
?
?0: (2.56)
29
In the following, we demonstrate that the perturbation formulations derived in sec-
tion 2.1 lead to the same conclusion.
2.2.1 Zero order solution
We flrst consider the zero order solution. The boundary value problem (2.14)-
(2.15) subject to a sphere (with radius r) can be rewritten in terms of the scalar
potential ? as follows:
r2?? = 0; (2.57)
?+ = ??; (2.58)
?(0)+ @?
+
@n = ?0
@??
@n ; (2.59)
and this potential satisfles the Sommerfeld radiation condition
limr!1pr
?@?
@r ?ik?
!
= 0: (2.60)
In this case, we are interested in such values of ?(0)+ that the sourceless flelds
may exist for the above boundary value problem. The general solutions of (2.55)
for regions V + and V ? are:
?+ (r; ;`) = X
l;m
AlmrlYlm ( ;`); (2.61)
??(r; ;`) = X
l;m
Blmr?l?1Ylm ( ;`); (2.62)
where Ylm ( ;`) are spherical harmonics, and Alm and Blm are unknown coe?cients.
Substituting (2.61)-(2.62) into boundary conditions (2.58)-(2.59), and by using or-
thogonality properties of spherical harmonics, we conclude that:
Alm = Blm; (2.63)
30
?(0)+ l = ??0 (l +1): (2.64)
Therefore,
?(0)+ = ?l +1l ?0: (2.65)
It is clear that when l = 1, the solutions of Laplace?s equation are three-
fold degenerate (m = 0;?1) and the associated electric flelds are spatially uniform
inside the sphere. These modes are the fundamental plasmon modes in spherical
nanoparticles and the corresponding resonance value of dielectric permittivity for
these modes is ?(0)+ = ?2?0. The frequency at which the dispersion relation assumes
this value is known as the Fr?ohlich frequency. To make the zero order solution
complete, it is necessary to know the scalar potential, surface charge density and
electrical flelds. Since the three modes are degenerate and the only difierence is
spatial orientation, in the following analysis, only the mode which has electric fleld
along the z-axis is considered. It is easy to show that (assuming A10 = 1):
?+ (r; ;`) =
s
3
4? cos ; (2.66)
??(r; ;`) =
s
3
4? cos ; (2.67)
 0 = @?
+
@n ?
@??
@n = 3
s
3
4? cos ; (2.68)
e+01 (r; ;`) = ?ar
s
3
4? cos ?a 
s
3
4? sin : (2.69)
In Cartesian coordinates, the electrical fleld for this mode can be expressed below
and it is clear that the electric fleld is uniform.
e+01x (x;y;z) = ?
s
3
4?ax (sin cos`cos ?cos cos`sin ) = 0; (2.70)
31
e+01y (x;y;z) = ?
s
3
4?ay (sin sin`cos ?cos sin`sin ) = 0; (2.71)
e+01z (x;y;z) = ?
s
3
4?az (cos cos +sin sin ) = ?
s
3
4?az: (2.72)
2.2.2 Second order solution
Since ?(0)+ = 2?0, the formula (2.54) for the second order correction of resonance
values of dielectric permittivity becomes:
?(2)+ = ?3
H
S nM ?~e
+
2 (M)?
+
0 (M)dSMH
S nM ?e
+
0 (M)?
+
0 (M)dSM
: (2.73)
Although the surface dipole density ?+0 is not explicitly computed in the zero order
solution, due to the fact that the resulting electric flelds are uniform for the mode
under consideration, ?+0 is proportional to  +0 up to a constant. Therefore,  +0 can
be used in the evaluation of the integrals in above formula. We flrst calculate the
denominator of formula (2.73):
Id =
I
S
nM ?e+0 (M) +0 (M)dSM: (2.74)
Substituting (2.68) and (2.69) into (2.74), after simple algebra, the analytical result
is obtained:
Id = ? 94?
Z 2?
0
d`
Z ?
0
cos2  r2 sin d = ?3: (2.75)
Then formula (2.73) becomes:
?(2)+ =
I
S
nM ?~e+2 (M)?+0 (M)dSM: (2.76)
Let?s denote this integral as In. The evaluation of integral In is somewhat lengthy,
nevertheless the main steps are specifled as follows. First, we need to express ~e+2 in
32
terms of e+0 . According to equation (2.50) and (2.69), it can be shown that:
~e+2 (M) = fi
I
S
brQM ?a`Q sin dSQ; (2.77)
where:
fi =
 
?+0
?0 ?1
?
8?
s
3
4? = ?
3
8?
s
3
4?: (2.78)
According to equation (2.77), formula (2.74) can be expressed as:
In = fi
I
S
arM ?
(I
S
rQ ?rM
rMQ ?(a`Q sin Q)dSQ
)
 +0 (M)dSM: (2.79)
Note for the unit sphere that
rM = arM = nM; rQ = arQ = nQ; (2.80)
which imply that formula (2.79) can be reduced to:
In = ?fi
I
S
a Q sin Q ?
(I
S
arM
rMQ 
+
0 (M)dSM
)
dSQ: (2.81)
The integral inside the bracket of (2.81) is a vector, the three components can be
analytically evaluated and the results are:
Ix = k1 cos Q sin Q cos`Q; (2.82)
Iy = k1 cos Q sin Q sin`Q; (2.83)
Iz = k1
 
cos2  Q + 43
?
; (2.84)
where
k1 = 65p3?: (2.85)
Substituting (2.82)-(2.84) into (2.81 ), we get:
In = ?fik1
I
S
?
cos2  Q sin2  Q ?sin2  Q
 
cos2  Q + 43
??
dSQ: (2.86)
33
The rest of the work is straightforward, with some caution in the evaluation, we
obtain:
?(2)+ = In = ?35: (2.87)
Then the resonance values of dielectric permittivity for spherical nanoparticles up
to second order corrections are found as
?+ = ?(0)+ +fl2?(2)+ = ?
 
2+ 35fl2
?
?0: (2.88)
This formula is identical to the analytical expression given by Mie theory (see for-
mula (2.54)). In conclusion, as far as the spherical nanoparticle is concerned, our
integral equation technique agrees with the classical Mie theory.
2.3 Numerical analysis of plasmon resonances in nanoparticles
2.3.1 Numerical technique for solution of integral equation
From the discussion in Section 2.1, it is clear that the key step of numerical
analysis of plasmon resonances in nanoparticles is to solve integral equation (2.21).
Now we proceed to the discussion of an e?cient numerical technique for the solutions
of this integral equation. To this end, let us partition S into N small pieces 4Sj
and rewrite integral equation (2.21) as follows:
 (Q) = ?2?
NX
j=1
Z
4Sj
 (M) rMQ ?nQr3
MQ
dSM: (2.89)
Now, we integrate (2.89) over 4Si:
Z
4Si
 (Q)dSQ = ?2?
NX
j=1
Z
4Sj
 (M)
"Z
4Si
rMQ ?nQ
rMQ dSQ
#
dSM; (i = 1;2;:::N):
(2.90)
34
By introducing the notation
!i (M) =
Z
4Si
rMQ ?nQ
r3MQ dSQ; (2.91)
the last formula can be presented as follows:
Z
4Si
 (Q)dSQ = ?2?
NX
j=1
Z
4Sj
 (M)!i (M)dSM: (2.92)
It is apparent that !i (M) is the solid angle which 4Si subtends at point M . By
introducing new variables
Xi =
Z
4Si
 (Q)dSQ; (2.93)
integrals in the right-hand side of (2.92) can be approximated as follows:
Z
4Sj
 (M)!i (M)dSM ? !i (Mj)Xj = !ijXj; (2.94)
whereMj issomemiddlepointofpartition4Sj. Itisapparent(onintuitivegrounds)
that approximation (2.94) is more accurate than direct discretization of integral in
(2.21), because solid angles !i (M) are smooth functions of M, while the kernel of
integral equation (2.21) is (weakly) singular. By substituting formulas (2.93) and
(2.94) into equation (2.92), we obtain:
Xi = ?2?
NX
j=1
!ijXj: (2.95)
Another advantage of discretization (2.95) is that the evaluation of singular integrals
in calculations of !ii can be completely avoided. Indeed, according to formulas (2.29)
and (2.91), we flnd:
!ii ? 2? ?
NX
i=1;i6=j
!ij: (2.96)
35
The numerical technique based on discretization (2.95) has been software imple-
mented and extensively tested by the author. It has proved to be remarkably accu-
rate, even for calculation of large eigenvalues.
2.3.2 Numerical results
The algorithm described above has been flrst tested for spherical nanoparticles
where exact analytical solutions are available (Mie theory[39]). The results of nu-
merical computations are presented in Table 2.1. It is apparent from this table that
numerical results are quite accurate even for appreciably high mode orders. The
surface charge distribution (eigenfunction) for the flrst 10 eigenmodes are plotted in
Figure 2.5. Next, the described algorithm has been tested for ellipsoidal nanoparti-
cles. The computational results are presented in Table 2.2 for the case of ellipsoid of
revolution with the main axis ratios 1:1:1.55. It is evident from this Table that the
computed eigenvalues compare quite well with the exact (theoretical) eigenvalues
?i = 11?2Ni for spatially uniform modes. It is also apparent that for those modes
P
i
1
?i is very close to 1 as it must be. The surface charge distribution (eigenfunction)
for the flrst 10 eigenmodes are plotted in Figure 2.6
Fig. 2.8 presents the computational results for the resonant free-space wave-
length as a function of separation distance between two gold spherical nanoparticles
located on a dielectric substrate with ? = 2:25?0 for difierent values of radius ratio
(see Fig. 2.7). It is clear from Fig. 2.8 that the separation between two spheres
can be efiectively used for tuning of plasmon resonances to desirable frequencies.
36
Fig. 2.8 is an example of our numerous computations performed for two and several
nanospheres of various radii and separation distances. We have computed the reso-
nance wavelength for two short gold cylinders of ellipsoidal cross-sections placed on
a dielectric substrate (see Fig. 2.9). Fig. 2.10 presents the computational results
for the resonance wavelength as the function of separation between cylinders for
difierent values of the axis ratio.
Table 2.3 presents the computational results for gold nano-rings placed on a
dielectric substrate. In this table, the computational results for resonance wave
lengths are compared with those found experimentally (see reference [72]). We
have also compared our numerical results with available experimental data for the
following two cases: a) a short gold cylinder of ellipsoidal cross-section (with long
axis 130nm, short axis 84nm, height 30nm) placed on a dielectric substrate and
b) a gold triangular nano-prism (with edge length 48nm, height 14nm). Table 2.4
presents the comparison between our computational results and experimental data
published in references [42, 73], respectively. The electrical fleld intensity of the flrst
two plasmon resonance modes for gold triangle prism is shown in Fig.2.11.
Figure 2.12 presents the comparison of computational and experimental results
[74] for the gold nanocube. For the case of the nanocube in water (Fig. 2.12(A)),
the two resonance peaks are observed in experiments at 432nm and 500nm, while
our computational results reveal resonances at 420nm and 478nm, respectively. For
the case of the nanocube on substrate (? = 2:25?0 ) (Fig. 2.12(B)), experiments
result in the two resonance peaks at 395nm and 457nm, while our computations
suggest resonances at 410nm and 443nm, respectively.
37
In all these computations, the experimentally measured dispersion relation for
gold [68] has been used. The meshes used in calculations are shown in Figure 2.13.
It is thus demonstrated that our computational results compare quite well with
known theoretical results and are in a reasonably good agreement with the available
experimental data.
Table 2.1: Eigenvalues for a single nanosphere.
Mode Computed Mie Mode Computed Mie
number eigenvalues theory number eigenvalues theory
1 2.999191 3 21 9.038890 9
2 2.999193 3 22 9.049236 9
3 2.999194 3 23 9.049301 9
4 4.980130 5 24 9.049312 9
5 4.980130 5 25 10.86267 11
6 4.980148 5 26 10.86268 11
7 5.022817 5 27 10.86291 11
8 5.022828 5 28 10.94108 11
9 6.927911 7 29 10.94200 11
10 6.981790 7 30 11.03838 11
11 6.981884 7 31 11.03839 11
12 6.981885 7 32 11.03846 11
13 7.027287 7 33 11.06465 11
14 7.027289 7 34 11.06497 11
15 7.027393 7 35 11.06500 11
16 8.915480 9 36 12.75603 13
17 8.915606 9 37 12.84736 13
18 8.915633 9 38 12.84819 13
19 8.979679 9 39 12.84824 13
20 8.979774 9 40 12.93150 13
38
Table 2.2: Eigenvalues for a single nanoellipsoid with aspect ratio 1:1:1.55.
Mode Computed Theoretical Mode Computed
number eigenvalues values number eigenvalues
1 1.976541 1.9723 6 4.590262
2 3.335450 7 4.829570
3 4.080630 4.0507 8 5.645674
4 4.080692 4.0570 9 5.645746
5 4.590155 10 6.314742
?4
?2
0
2
4
?4
?2
0
2
4
?5
0
5
?2
0
2
4
6
?5
0
5
?5
0
5
?5
0
5
?5
0
5
?5
0
5
?5
0
5
Figure 2.5: Surface charge distribution (eigenfunction) for the flrst 10 eigenmodes
of single nanosphere.
39
?5
0
5
?2
0
2
4
6
?2
0
2
?2
0
2
?4
?2
0
2
4
?4
?2
0
2
4
?5
0
5
?5
0
5
?5
0
5
?8
?6
?4
?2
0
2
Figure 2.6: Surface charge distribution (eigenfunction) for the flrst 10 eigenmodes
of single nano-ellipsoidal.
40
Table 2.3: Comparison with experimental results [72] for gold nanoring.
Ring1 Ring2 Ring3
Outer radius of the ring (nm) 60 60 60
Height of the ring (nm) 40 40 40
Thickness of ring wall (nm) 14(2) 10(2) 9(2)
Experimental resonance (nm) 1000 1180 1350
Computational resonance (nm) 940 1102 1214
Table 2.4: Comparison with experimental results [42, 73] for gold ellipsoidal cylinder
and triangular prism.
Resonance cylinder prism
wavelength (nm) (nm)
Computational results 622 653
Experimental results 645 690
41
DR1 R2
substrate
Figure 2.7: Two nanospheres on substrate.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2490
500
510
520
530
540
550
560
570
580
590
Resonance wavelength (nm)
Separation/Radius
1 
2 
3 
4 
5 
             R1:R2curve1 ? 1.0:1.0  
curve2 ? 1.0:0.9  curve3 ? 1.0:0.8  
curve4 ? 1.0:0.7  curve5 ? 1.0:0.6  
Figure 2.8: Resonance wavelength for two nanospheres on substrate.
42
DxD xD
Figure 2.9: Two ellipsoidal cylinders on substrate.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1550
600
650
700
750
800
Resonance wavelength (nm)
Separation/Short axis length
          Dx:Dy       curve1 ? 2.0:1.0      
curve2 ? 1.8:1.0      curve3 ? 1.6:1.0      
curve4 ? 1.4:1.0      curve5 ? 1.2:1.0      
curve6 ? 1.0:1.0      
1 
2 
3 
4 
5 
6 
Figure 2.10: Resonance wavelength for two ellipsoidal cylinders on substrate.
43
Figure 2.11: The plasmon mode light intensity of flrst two resonance modes for
triangular nanoprism, colorbars indicate the ratio of optical near fleld intensity nor-
malized to incident flelds.
44
Figure 2.12: Comparison of the plasmon resonance wavelength of (A) nanocube
ensemble in water and (B) single nanocube on a glass substrate. The down arrows
indicate the experimental results (432 nm and 452 nm in (A), 395 nm and 457 nm
for (B)) while the up arrows show the computational results (420 nm and 478 nm
in (A), 410 nm and 443 nm in (B)).
45
(a) sphere (b) cube
(c) ring (d) triprism
Figure 2.13: Examples of surface mesh used in computations for difierent geometry
of three dimensional nanoparticles.
46
2.4 Plasmon resonances in shell structures
The tunability of plasmon resonant frequencies in nanoparticles is important
for many applications such as nanophotonics and biosensors. In section 2.3, we
have seen that the plasmon resonant frequencies can be tuned by using a two-sphere
system and adjusting the coupling between the two particles. The tunable frequency
range of this conflguration is limited, however, and it usually experiences more losses
than single nanoparticle. The wide range tunability of plasmon resonant frequencies
of single nanoparticle is clearly desired. It has been observed [27, 28] that this
feature can be achieved by using metallic nanoshells and controlling the resonances
frequencies via the shell thickness. In this section, the boundary integral equation
technique for direct calculations of resonance frequencies of metallic nanoshells is
presented. This technique is a modiflcation of the method developed in previous
sections for solid nanoparticles. The wide range tunability of plasmon resonant
frequencies is illustrated by numerical examples. In addition, analytical studies are
performed for spherical and ellipsoidal shells, and a phenomenon of twin spectrum
is discussed.
2.4.1 A generalized eigenvalue problem
To start the discussion, consider a metallic nanoshell with dielectric core shown
in Figure 2.14(a). The plasmon resonances in metallic nanoshell structures can be
formulated as the following boundary value problem: flnd such values of ?+ for which
47
( )? ?+

?
?








?
?
?









?
?
 
Figure 2.14: (a) Schematic of a nanoshell with dielectric core; (b) the equivalent
problem for plasmon resonance in nanoshells.
there exist non-zero solutions to the difierential equations:
r2?1 = 0 in V1; (2.97)
r2?2 = 0 in V2; (2.98)
r2?3 = 0 in V3; (2.99)
subject to the boundary conditions:
?1@?
+
1
@n = ?+
@?+2
@n ; ?
+ = ??; on S1; (2.100)
?+@?
+
2
@n = ?0
@?+3
@n ; ?
+ = ??; on S2; (2.101)
?3(1) = 0: (2.102)
Here: ?1 and ?+(!) are dielectric permittivities for the core and metallic nanoshell,
respectively.
The electric potential of this sourceless fleld can be represented as an electric
potential of single layers of electric charges  1 and  2 distributed over S1 and S2,
48
respectively:
?(Q) = 14??
0
"I
S1
 1(M)
rMQ dSM +
I
S2
 2(M)
rMQ dSM
#
: (2.103)
In other words, single layers of electric charges  1(on S2) and  2(on S2) may create
the same electric fleld in homogeneous space with ?0 as the resonant source-free
electric fleld that may exist in the presence of the nanoshell and dielectric core
(see Figure 2.14(b)). It is apparent that the potential given by the formula (2.103)
satisfles the Laplace equations (2.97)-(2.99) and continuity conditions on S1 and
S2. By using formula (2.20), it can be concluded that the boundary conditions
(2.100)-(2.101) for normal derivatives will be satisfled if  1 and  2 are solutions of
the following coupled homogeneous boundary integral equations:
?1
?
2? 1(Q)+HS1  1(M)rMQ?nQr3
MQ
dSM +HS2  2(M)rMQ?nQr3
MQ
dSM
?
= ?+
?
?2? 1(Q)+HS1  1(M)rMQ?nQr3
MQ
dSM +HS2  2(M)rMQ?nQr3
MQ
dSM
?
;
(2.104)
?0
?
?2? 2(Q)+HS1  1(M)rMQ?nQr3
MQ
dSM +HS2  2(M)rMQ?nQr3
MQ
dSM
?
= ?+
?
2? 2(Q)+HS1  1(M)rMQ?nQr3
MQ
dSM +HS2  2(M)rMQ?nQr3
MQ
dSM
?
:
(2.105)
Thus, source-free electric flelds may exist only for such negative values of permit-
tivity ?+(!) that the coupled homogeneous integral equations (2.104)-(2.105) have
non-zero solutions. By discretizing integral operators in the above integral equa-
tions, we arrive at the following generalized eigenvalue problem for ?+:
2
66
64
?1(K11 +2?I) ?1K12
?1K12 ?0(K22 ?2?I)
3
77
75
2
66
64
?!X
1
?!X
2
3
77
75 = ?+
2
66
64
K12 ?2?I) K12
K12 K22 +2?I
3
77
75
2
66
64
?!X
1
?!X
2
3
77
75;
(2.106)
49
where matrices Kij are discretized versions of the corresponding integral operators
in (2.104)-(2.105), ?!X2 and ?!X2 are discretized versions of  1(M) and  2(M), respec-
tively, while I is the identity matrix. The solution of the generalized eigenvalue
problem yields the resonance values of ?+, and then the appropriate dispersion rela-
tion ?+(!) can be employed to flnd the resonance frequencies for metallic nanoshells.
Since the integral operators in (2.104) and (2.105) are compact, the spectrum is dis-
crete. Furthermore it can be shown that all the eigenvalues associated with the
generalized eigenvalue problem (2.97)-(2.101) are real and negative.
It can be observed that the mathematical structure of boundary integral equa-
tion (2.104)-(2.105) is invariant with respect to the simultaneous and identical scal-
ing of S1 and S2, i.e., the scaling of the dimensions of the shell. This results in
the unique property of plasmon resonances: resonance frequencies depend on shell
shapes but they are scale invariant with respect to shell dimensions, provided that
these dimensions are appreciably smaller than the free-space wavelength.
Table 2.5: Comparison of the resonance wavelengths for Au spherical nanoshells
with silicon cores.
Shell1 Shell2
Inner/outer radius 60/80 nm 55/75nm
Experimental results [27] 752 nm 717 nm
Computational results 795 nm 764 nm
The above integral equation based generalized eigenvalue problem has been
solved by using the discretization technique described in Section 2.3. Table 2.5
presents the results for the plasmon resonances of gold spherical nanoshells with
50
silicon cores computed by using our technique and experimentally measured in [27].
The dispersion relation for gold published in [68] has been used in calculations. This
table reveals a fairly good agreement between the experimental and computational
results.
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3500
600
700
800
900
1000
1100
1200
1300
outer radius / inner radius
Resonance wavelength
Figure2.15: PlasmonresonantwavelengthsforAusphericalnanoshellswithdifierent
shell thickness. The core of these shells is silicon.
The resonance wavelength of gold spherical nanoshells (with silicon core) has
been also computed for various ratios of inner and outer radius. The numerical
results are illustrated in Figure 2.15. This flgure clearly reveals that the resonance
wavelength of nanoshell structures can be tuned within a wide range (from visible
to near infrared) by adjusting the thickness of the shell.
51
2.4.2 Analytical results for spherical and ellipsoidal shells
For some simple geometries, it is possible to pursue analytical solutions for
the boundary value problem (2.97)-(2.101). In the following, we demonstrate these
solutions for spherical and ellipsoidal nanoshells and discuss the \twin" spectra
phenomena associated with these shell structures. We also compare these analytical
results with numerical computations of the generalized eigenvalue problem.
For spherical shells, the solution of the boundary value problem (2.97)-(2.101)
can be represented in terms of the spherical harmonics:
?1(r; ;`) = ArlYlm( ;`); (2.107)
?2(r; ;`) = [Brl +Cr?(l+1)]Ylm( ;`); (2.108)
?3(r; ;`) = Dr?(l+1)Ylm( ;`); (2.109)
where A;B;C and D are unknown coe?cients. By substituting (2.107)-(2.109) into
boundary conditions (2.100)-(2.101), the linear homogeneous equations for the four
unknown coe?cients are obtained. The occurrence of plasmon resonances implies
the existence of non-zero solutions to these equations. The existence of non-zero
solution, in turn, requires that the determinant of the above homogeneous equations
is equal to zero. In this way, the following quadratic equation for ?+ can be derived:
?2+l + (?1 +?0)l(l +1)+?1l
2(r2=r1)3 +?0(l +1)2(r2=r1)3
l(l +1)[(r2=r1)3 ?1] ?+l +?1?0 = 0; (2.110)
where r1 and r2 are the inner radius and outer radius of the spherical nanoshell,
respectively. It turns out that for the spherical nanoshell structure an interesting
phenomenon of twin spectra occurs. The essence of this phenomenon is that for any
52
Table 2.6: Comparison of the resonance values of ?+ for spherical nanoshells, ?1 =
5?0.
Shell1 Shell2
r2=r1 1.5:1 20:1
l = 1 (Theoretical) -8.32,-0.60 -2.50,-1.99
l = 1 (Computational) -8.35,-0.60 -2.51,-1.99
l = 2 (Theoretical) -8.79,-0.58 -3.32,-1.51
l = 2 (Computational) -8.83,-0.57 -3.34,-1.50
l = 3 (Theoretical) -9.15,-0.56 -3.71,-1.34
l = 3 (Computational) -9.21,-0.55 -3.75,-1.33
resonant mode associated with index l there are two resonant values ?(1)+l and ?(2)+l
related by the formula:
?(1)+l = ?1?0?(2)
+l
: (2.111)
Indeed, formula (2.111) is the immediate consequence of equation (2.110). Table
2.6 presents the results for plasmon resonances in spherical nanoshells computed by
using the integral equation based numerical technique and the analytical formula
(2.110). This table illustrates an excellent agreement between the theoretical and
computational results. It is clear that the twin values of ?+l satisfy (2.111).
In the case of ellipsoidal nanoshells, spatially uniform resonance modes can be
easily found analytically. For these modes the relation between the polarization P
of the shell and applied uniform fleld E0 can be expressed in the form:
2
66
66
66
66
4
fi+x 0 0
0 fi+y 0
0 0 fi+y
3
77
77
77
77
5
P = E0: (2.112)
53
By using ellipsoidal coordinates [39], the following formulas for fi+x , fi+y and fi+z can
be derived:
fi+i =
h
?+ +(?1 ??+)
?
N(1)i ?fN(2)1
?ih
?0 +(?+ ??0)N(2)i
i
+fN(2)i ?+(?1 ??+)
?
n
(?+ ??0)
h
?+ +(?1 ??+)
?
N(1)i ?fN(2)1
?i
+f?+(?1 ??+)
o ;
(2.113)
where i = x;y;z;
? = 4?3 a2b2c2; f = a1b1c1a
2b2c2
; (2.114)
a1;b1;c1 are semi-axes of the inner ellipsoid, a2;b2;c2 are semi-axes of the outer
ellipsoid, while depolarization factors Ni are given by the formulas:
Nx = abc2
Z 1
0
d?
(?+a2)
q
(?+a2)(?+b2)(?+c2)
; (2.115)
Nx = abc2
Z 1
0
d?
(?+b2)
q
(?+a2)(?+b2)(?+c2)
; (2.116)
Nx = abc2
Z 1
0
d?
(?+c2)
q
(?+a2)(?+b2)(?+c2)
; (2.117)
and Nx+Ny +Nz = 1. In the case of axially symmetric ellipsoids (b = c), the above
formulas are reduced to:
Nx = 1?fl
2
e2
" 1
2flln
?1+fl
1?fl
!
?1
#
; Ny = Ny = 1?Nx2 ; (2.118)
fl2 = 1? b
2
a2: (2.119)
The existence of resonances implies that the matrix in (2.112) is singular. According
to (2.113), this means that the resonance values of ?+(!) can be computed by solving
the following quadratic equation:
?
1?f ?N(1)i +fN(2)i
?
N(2)i ?2+ +
?
N(1)i ?fN(2)i
??
1?N(2)i ?1?0
?
+
h?
1?N(1)i +fN(2)i
??
1?N(2)i
?
+
?
f +N(1)i ?fN(2)i
?
N(2)?1i
i
?+ = 0:
(2.120)
54
If the shell has self-similar boundaries S1 and S2, the inner ellipsoid and the outer
ellipsoid have the same depolarization factors. Then, formula (2.120) is reduced to:
?2+ + [(1?Ni +fNi)(1?Ni)?0 +(f +Ni ?fNi)Ni?1](1?f)(1?N
i)Ni
?+ +?1?0 = 0: (2.121)
Equation (2.121) has the same structure as formula (2.110). This suggests that
the phenomenon of twin spectra occurs for ellipsoidal shells as well, albeit only for
uniform modes. The two resonant values ?(1)+ and ?(1)+ are related by the formula:
?(1)+ = ?1?0?(2)
+
: (2.122)
Table 2.7 presents the results for plasmon resonances in ellipsoidal nanoshells com-
puted by using the integral equation based numerical technique and the analytical
formula (2.121). This table illustrates a very good agreement between the theoretical
and computational results.
Table 2.7: Comparison of the resonance values of ?+ for ellipsoidal nanoshells, ?1 =
15:2?0.
Shell1 Shell2
a : b : c 1.414:1:1 1.414:1:1
Shell ratio 4:3 7:5
?+(Theoretical) -25.00,-0.61 -21.04,-0.71
?+(Computational) -25.83,-0.59 -22.08,-0.69
In conclusion, the analysis of plasmon resonances in metallic nanoshells with
dielectric cores can be reduced to an generalized eigenvalue problem. For spherical
and ellipsoidal nanoshells, simple analytical expressions are obtained for the direct
calculation of plasmon resonant frequencies. The wide range tunability of plasmon
55
resonant frequencies via adjusting the thickness of shells has been illustrated by
numerical examples. Furthermore, it can be seen from formulas (2.106), (2.110) and
(2.120), the dielectric permittivity of the core also plays a role in the determination
of the plasmon resonant frequencies, which adds additional capability to tune the
resonant frequency.
56
Chapter 3
Fast Multipole Method for the Solution of Eigenvalue
Problem
In this chapter, the fast multipole method (FMM) is introduced for the solution of
large scale eigenvalue problem. The applicability of FMM and its implementation in
the iterative algorithm are presented in detail. Numerical results which illustrate the
e?ciency of FMM compared with traditional method for the numerical analysis of
plasmon resonances are presented. In addition, a semi-analytical method based on
multipole expansions is developed for the analysis of plasmon resonances in spherical
nanoparticles.
57
3.1 Why fast multipole method is needed?
In Chapter 2, we have demonstrated that the plasmon resonance frequencies
can be directly obtained by solving the following homogeneous boundary integral
equations
 (Q) = ?2?
I
S
 (M) rMQ ?nQr3
MQ
dSM; (3.1)
where ? is given by (2.22). By using the discretization technique described in section
2.3, this eigenvalue problem can be written in the matrix form as
Xi = ?2?
NX
j=1
!ijXj; (3.2)
where Xi and !ij are deflned in (2.89) and (2.91), respectively. The solutions of the
eigenvalue problem (3.2) is usually obtained by using iterative techniques, which
require matrix-vector product B ? Avk computations. It is apparent that when
matrix A is fully populated, the multiplication Avk is accomplished in O(N2) op-
erations. For large N, the computation cost can be prohibitive. In order to solve
this problem, we have to seek a more advanced algorithm. It is turns out that fast
multipole method (FMM) is ideally suitable for this situation and can appreciably
save the cost of CPU time and memory.
The fast multipole method is regarded as one of the greatest algorithms of
the 20th century [65, 75, 76]. This algorithm allows the product of particular dense
matrices with a vector to be evaluated approximately (to a specifled precision) in
O(NlogN) operations, when direct multiplication requires O(N2) operations (when
N is large, the save is considerable). The FMM is originally proposed by Rokhlin
[77] as a fast scheme for accelerating the numerical solution of the Laplace equation
58
in two dimensions. It was further improved by Greengard and Rokhlin [64, 79] when
it was applied to particle simulations. Since then, FMM has wide applications in
astrophysics, molecular dynamics, computational  uid dynamics and radar scatter-
ing, etc. For example, Mayergoyz et al applied FMM in the nonlinear magnetosta-
tic calculations [78] recently. Furthermore, Gumerov and Duraiswami published a
monograph [65] in 2005 which reveals the most recent advances in FMM. Its cen-
tral strategy is the hierarchical decomposition of the data-space in the form of a
quadtree (or octtree for the 3-dimensional case). This hierarchical decomposition is
used to cluster particles in the computational domain at various spatial lengths and
compute interactions with other clusters that are su?ciently far away by means of
series expansions (see Figures 3.1 and 3.2).
Figure 3.1: Space partitioning. Figure 3.2: Data structure.
FMM does not work for any matrix-vector multiplication, instead, it only
works for cases that satisfy the following conditions: (1) two sets of points (for ex-
ample, source points and evaluation points) are involved in a vector space; (2) the
59
dense matrix is generated by potential (kernel) functions; (3) these functions can be
factorized and have local and multipole expansions; (4) the local expansion coe?-
cients can be RjR-translated (i.e., local-to-local translation) and the multipole ex-
pansion coe?cients can be SjS-translated (i.e., multipole-to-multipole translation);
(5) The multipole expansion coe?cients can be SjR-translated (i.e., multipole-to-
local translation). Since the matrix elements of (3.2) are derived from a 1=r-type
kernel, it can be shown that all the above conditions above are met.
3.2 Factorization and translations of the kernel function
In this section, we examine the applicability of FMM to our 1=r-type kernel
function from the mathematic point of view. The key components such as factor-
ization, expansion and translation are described below [65, 80].
3.2.1 Factorization of the kernel function
The kernel function of integral equation (3.1) is:
!ij =
Z
?Si
rMjQi ?nQi
r3MjQi dSi: (3.3)
If the patch is small enough, we can assume the distribution of the physical quantity
is uniform on the patch, the integral in the last equation becomes a multiplication
!ij = rMjQi ?nQir3
MjQi
?Si; (3.4)
and it is indeed the gradient of 1=rMQ:
!ij = r( 1r
MjQi
)?nQi?Si; (3.5)
60
where
rMjQi =
flfl
flrMj ?rQi
flfl
fl: (3.6)








Figure 3.3: Multipole expansion.
Now we start from the factorization of kernel 1r. The expansion in the three-
dimensional case can be obtained according to the additional theorem of spherical
harmonics:
1
jx?x0j =
1X
l=1
lX
m=?l
4?
2l +1
rl<
rl+1> Y
?
lm( 
0;`0)Ylm( ;`); (3.7)
where
Y ?lm( ;`) = (?1)mYl;?m( ;`); (3.8)
the x0(r0; 0;`0) and x(r; ;`) indicate the source point and evaluation point respec-
tively. Supposing we select the expansion center at x?(r1; 1;`1), we can introduce
the multipole expansion (S-expansion) and local expansion (R-expansion). Let x
be the evaluation point and x0 be the source point. If jx?x?j > jx0 ?x?j, the
multipole expansion is deflned as the following (see Figure 3.3):
1
j(x?x?)?(x0 ?x?)j =
1X
l=1
lX
m=?l
Clm(x0?x?)Slm(x?x?); (3.9)
61
where Clm(x0?x?) = 4?2l+1 jx0 ?x?jl Y ?lm(x0 ?x?); Slm(x?x?) = Ylm(x?x?)jx?x
?jl+1
.








Figure 3.4: Local expansion.
For the case jx0 ?x?j > jx?x?j, the R-expansion is deflned as following:
1
j(x?x?)?(x0 ?x?)j =
1X
l=1
lX
m=?l
Blm(x0?x?)Rlm(x?x?); (3.10)
where Blm(x0?x?) = 4?2l+1 Y?lm(x0?x?)jx0?x
?jl+1
; Rlm(x?x?) = jx?x?jl Ylm(x?x?).
3.2.2 Translation in 3D
In FMM algorithm, the multipole (S-) expansion and local (R-) expansion
of a potential need to be translated to the same type of expansion at difierent
expansion centers or a difierent type of expansion at difierent expansion centers.
In the following, we discuss three types of translation: the \multipole-to-multipole"
translation, the \local-to-local" translation and the \multipole-to-local" translation.
\Multipole-to-multipole" translation. Consider the multipole expan-
sion in formula (3.9) with expansion center x?(r1; 1;`1), this expansion can be
62
Figure 3.5: Illustration of local-to-local (RjR), multipole-to-local (SjR), and
multipole-to-multipole (SjS) translations from expansion center r1 to expansion
center r2 (translation vector t = r2 ? r1). The star shows location of a source.
Expansion about center r1 is valid inside the lighter and darker gray area, while
expansion about center r2 is valid only in the darker gray area.
re-expanded with the center at x0?(r2; 2;`2) and expressed as:
1
jx?x0j =
1X
j=0
jX
k=?j
Mjk(x0 ?x0?)Sjk(x?x0?); (3.11)
where
Mjk =
1X
n=0
nX
m=?n
Ok?mj?n Jk?mm Amn Ak?mj?n
Akj jx
0 ?x0
?j
n Y ?m
n (x
0 ?x0
?); (3.12)
Amn = (?1)
n
q
(n?m)!(n+m)!
; (3.13)
Jm0m =
8
>>><
>>>:
(?1)min(jm0j;jmj); m?m0 < 0;
1; otherwise:
(3.14)
\Local-to-local" translation. Consider the local expansion in formula (3.10) with
63
expansion center x?(r1; 1;`1), this expansion can be re-expanded with the center
at x0?(r2; 2;`2) and expressed as:
1
jx?x0j =
1X
j=0
jX
k=?j
Ljk(x0 ?x0?)Rjk(x0 ?x0?); (3.15)
where
Ljk =
1X
n=0
nX
m=?n
Omn Jmn?j;m?kAkjAm?kn?j
Amn jx
0 ?x0
?j
n?j Y m?k
n?j (x
0 ?x0
?); (3.16)
Jm0n;m =
8>
>>>>
>>><
>>>
>>>>
>:
(?1)n(?1)m; m?m0 < 0;
(?1)n(?1)m0?m; m?m0 < 0andjm0j < jmj;
(?1)n; otherwise;
(3.17)
and Amn deflned in (3.13).
\Multipole-to-local" translation. Consider the multipole expansion in formula
(3.9) with expansion center x?(r1; 1;`1), this expansion can be re-expanded to a
local expansion with the center at x0?(r2; 2;`2) and expressed as:
1
jx?x0j =
pX
j=0
jX
k=?j
Ljk(x0 ?x0?)Rjk(x0 ?x0?); (3.18)
where
Ljk =
pX
n=0
nX
m=?n
Omn (?1)n+kAmn Akj
Am?kn+j jx
0 ?x0
?j
?(n+j+1) Y m?k
n+j (x
0 ?x0
?); (3.19)
with Amn deflned in (3.16).
The translations presented above are inflnite in dimension and the translations
are exact operations. However, in the actual FMM algorithm, a truncated matrix
of size p2 ? p2 (for 3D case) is used instead of the actual transformation matrix,
64
thus introducing an error in the computation. This truncation is desirable because
it reduces the computational complexity, and it can be shown that [64, 65], for
convergent series, the error introduced in this step is bounded and can be reduced
to an arbitrarily small number by increasing the truncation number p. It is also
possible to determine what truncation number p to use in order to achieve a desired
level of accuracy [65, 80].
3.2.3 Difierentiation of the 1r kernel in 3D
Now we consider the gradient of kernel 1=rMQ since the real computation relies
on it. After factorization, the gradient is carried out on each item of the series. So
we need to know how to calculate the gradient of the spherical harmonics. Since
the data structure and the space discretization requires the Cartesian coordinates,
we also need to convert the difierentiation of spherical harmonics from spherical
coordinates to the Cartesian coordinates. At the same time, we will employ the
recursive relation of Legendre functions to get the expression. The derivation is
straightforward but somewhat lengthy, in the following, we will only give the results
[65].
In spherical coordinates, the gradient is deflned as
rx
? 1
jyj ?xij
!
=
"
ar @@r
x
+a @r
x@ 
+a` @r
x sin x`
#
; (3.20)
1X
n=0
nX
m=?n
4?
2n+1
rn<
rn+1> Y
?
nm( xi;`xi)Ynm( yj;`yj): (3.21)
The conversions between two coordinate systems are given by
x = rsin cos`; y = rsin sin`; z = rcos ; (3.22)
65
? = cos ; (3.23)
Dx = 12
?
Dxy +Dxy
?
(3.24)
Dy = i2
?
Dxy ?Dxy
?
; (3.25)
Dz = @@z = ? @@r + 1??
2
r
@
@?; (3.26)
where
Dxy = @@x +i @@y = e
i`
rp1??2
"
(1??2)(r @@r ?? @@?)+i @@`
#
; (3.27)
Dxy = @@x ?i @@y = e
?i`
rp1??2
"
(1??2)(r @@r ?? @@?)?i @@`
#
: (3.28)
The recursive relations of the associated Legendre functions are
?Pmn = 12n+1
h
(n+m)Pmn?1 +(n?m+1)Pmn+1
i
; (3.29)
(1??2) dd?Pmn = 12n+1
h
(n+1)(n+m)Pmn?1 ?n(n?m+1)Pmn+1
i
: (3.30)
Finally, we can get:
Dz (Rmn ) = amn Rmn?1; (3.31)
Dz (Smn ) = ?amn+1Smn+1; (3.32)
where:
amn =
s
2n+1
2n?1(n
2 ?m2); (3.33)
Rmn = rnY mn ( ;`); (3.34)
Smn = r?n?1Y mn ( ;`); (3.35)
Y mn ( ;`) = (?1)m
vu
ut2n+1
4?
(n?jmj)!
(n+jmj)!P
jmj
n (?)e
im`; (3.36)
66
Dxy(Rmn ) = bmn Rm+1n?1 ; (3.37)
Dxy(Smn ) = c?m?1n+1 Sm+1n+1 ; (3.38)
Dxy(Rmn ) = b?mn Rm?1n?1 ; (3.39)
Dxy(Smn ) = cm?1n+1 Sm?1n+1 ; (3.40)
where
bmn =
8
>>>>
>>>>
<
>>>>
>>>
>:
?
r
(2n+1)(n?m?1)(n?m)
(2n?1) ; 0 ? m ? n;
r
(2n+1)(n?m?1)(n?m)
(2n?1) ; ?n ? m ? 0;
0; jmj > n
(3.41)
cmn =
8>
>>>>
>>><
>>>
>>>>
>:
r
(2n?1)(n?m?1)(n?m)
(2n+1) ; 0 ? m ? n;
?
r
(2n?1)(n?m?1)(n?m)
(2n+1) ; ?n ? m ? 0;
0; jmj > n:
(3.42)
Dx(Rmn ) = 12[bmn Rm+1n?1 +b?mn Rm?1n?1 ]; (3.43)
Dx(Smn ) = 12[c?m?1n+1 Sm+1n+1 +cm?1n+1 Sm?1n+1 ]; (3.44)
Dx(Rmn ) = i2[bmn Rm+1n?1 ?b?mn Rm?1n?1 ]; (3.45)
Dx(Smn ) = i2[c?m?1n+1 Sm+1n+1 ?cm?1n+1 Sm?1n+1 ]: (3.46)
So, the basic conclusion is that those multipole or local expansion coe?cients can
be expressed through recursive relations. This is very helpful for software imple-
mentation.
67
3.3 FMM algorithm and implementation
3.3.1 FMM algorithm
Although the detailed description of FMM algorithm can be found in many
references [65, 67, 79], in order to make this dissertation self-complete, the FMM
algorithm is stated here in a short format. Before starting, we introduce some
notation. 'l;i is the p?term multipole expansion about the center of the box i at
level l, describing the potential fleld outside box i?s neighbors due to all particles
contained inside the box i. ?l;i is the p?term local expansion about the center of
the box i at level l, describing the potential fleld induced by all particles outside box
i?s neighbors. e?l;i is the p?term local expansion about the center of the box i at
level l, describing the potential fleld induced by all particles outside the neighbor?s
of box i?s parent.
Initialization. Choose precision to be desired ? and the number of levels n
for the hierarchy data structure. Set up the hierarchical data structure and sort all
points into the boxes at the flnest level.
Upward pass
Step 1
For each box i at the flnest level n, form the p-term multipole expansion
'n;i, representing the potential fleld induced by all particles in the box. record the
coe?cients of each expansion for all boxes.
Step 2
For levels l = n ? 1;n ? 2;:::;2, for each box i at level l, use multipole-to-
68
multipole translation to shift the centers of the multipole expansions of its children
boxes to the center of box j and merge them to form 'l;j which represents the
potential fleld induced by all particles in box j, or all particles in its children boxes.
See Fig.3.6.
Figure 3.6: Step 2 of FMM algorithm: Multipole to multipole translation.
Downward pass
Step 3
For levels n = 2;3;:::;n, for each box j at level l, use local-to-local translation
(See Fig.3.7). to shift the center of local expansion ?l?1 of j?s parent box to the
center of box j to form e?l;i. Use multipole-to-local translation (See Figs. 3.8 and
3.9) to convert the multipole expansion of all boxes that is in the interaction list of
box j to a local expansion about the center of box j, and add all these to e?l;i to
form ?l;i, which represents the potential fleld induced by all particles outside box
69
j?s neighbors.
Figure 3.7: Step 3 of FMM algorithm: Local to Local translation.
Step 4
For all the boxes j at the flnest level, for each particles at box j, evaluate 'n;j
at the particles position.
Step 5
For all the boxes j at the flnest level, for each particles at box j, evaluate the
interactions with particles in j?s neighbor boxes directly. The  owchart of above
described FMM algorithm can be seen in Figure 3.10.
70
Figure 3.8: Step 3 of FMM algorithm: Multipole to local translation.
Figure 3.9: Step 3 of FMM algorithm: Multipole to local translation.
71
Figure 3.10: A  ow chart of the standard FMM.
3.3.2 Date structure for FMM implementation
A key issue in the implementation of fast multipole methods is to establish
a well-organized data structure. In order to exploit the observation that the efiect
of a cluster of particles at a certain distance may be approximated by a flnite sum
of series expansions using the equations described in previous sections, we need to
organize the particles in a hierarchy of clusters. This hierarchy of clusters allows one
to e?ciently determine when the approximation is valid and therefore correlates to
the overall performance of the FMM algorithm.
Since the regions of expansion validity are specifled in terms of Euclidean
distance, subdivision of space into d-dimensional cubes is convenient for range eval-
uation between points, which is usually called 2d-tree data structure. In the 3D
case, the hierarchy tree structure can be seen in Figure 3.11. The main technique
for working with 2d-trees (and k-d trees) is the bit-interleaving technique [81, 82].
This technique enables O(1), or constant, algorithms for parent and sibling search
and O(logN) algorithms for neighbor and children search.
72
Figure 3.11: Hierarchical cells at three levels in three dimensions.
To e?ciently determine the \Parent" box and \Children" box, we can index
the boxes in the following way. Each box in the tree can be identifled by assigning
it a unique index among the 2d children of its parent box, and by knowing the index
of its parent in the set of its grandparent?s children. Then the index of a box can
be written as the string
String(n;l) = (N1;N2;:::Nl); Nj = 0;:::2d ?1; j = 1;:::;l; (3.47)
where l is the level at which the indexed box is located and Nj is the index of a box
at level j containing that box. Now the indexing string can be converted to a single
number as follows:
n = (2d)l?1 ?N1 +(2d)l?2 ?N2 +2d ?Nl?1 +Nl: (3.48)
Note that this index depends on the level l at which the box is considered and unless
this information is included, difierent boxes could be described by the same index.
The unique index of any box can be represented by the pair:
UniversalIndex = (n;l): (3.49)
73
If the indexing at each level is performed in a consistent way, then we call such a
indexing scheme \hierarchical". A consistent hierarchical scheme has the following
desirable properties.
1. Determining the Parent: Consider a box at level l of the 2d-tree, whose
index is given by equation (3.48). It can be shown that the parent of this box is
Parent(n;l) = (Parent(n);l?1): (3.50)
2. Determining the Children: For the universal numbering system 3.49,to get
the indices of all 2d children of a box represented by the string 3.47, the operation of
flnding the children is simply the calculation of the children numbers and assigning
their level to l +1:
ChildrenAll(n;l) = (ChildrenAll(n);l +1): (3.51)
The use of 2d-trees makes obtaining parent and children indices very convenient.
Indeed, the above operations are nothing but shift operations in the bit represen-
tation of n. Performing a right bit-shift operation on n by d-bits one can obtain
the index of the parent. One can list all indices of the children boxes of n by a left
bit-shift operation on n by d-bits and adding all possible combinations of d-bits.
In addition to flnd the \parent" and \children" for a box, the FMM also needs
a way to determine neighbors and the box center for a given index (n;l) and the
box index to which a given spatial point x belongs. To do this, a spatial ordering
in d-dimensional space should be introduced. Below, such an ordering and O(1)
algorithms for these operations are introduced.
74
We flrst scale the computation domain (the largest box) to the unit cube.
Any point x in the original three-dimensional space now can be found given ?x 2
[0;1] ? [0;1] ? [0;1]. Since all the boxes are already ordered by their indices at a
given level l from 0 to 2l?1, there is a straightforward correspondence between box
indices and coordinates of points. For this reason, we have the following convenient
operations.
3. Finding the index of the box containing a given point. Consider the relation
between the coordinate of a point and the index of the box where the point is
located. We note that the size of a box at each level is 1 placed at the position equal
to the level number after the decimal in its binary record, and it can be shown that
the operation is simply done by:
(n;l) = [23l ? ?x]: (3.52)
where [] means integral part.
4. Finding the center of a given box. This operation can be done through
?xm(n;l) = 2?l ?(nl +2?1); m = 1;2;3: (3.53)
5. Finding neighbors. This operation to the k-neighbor can be done through
Neighborm =
8>
>>>
>>>>
<
>>>>
>>>>
:
n?m;nm;n+m; nm 6= 0;2l ?1;
nm;n+m; nm = 0;
n?m;nm; nm = 2l ?1:
(3.54)
where n+m = nm +k, n?m = nm ?k, m = 1;2;3.
The data structure and the indexing and ordering technique described above
75
have been implemented in code and it is proven to be very e?cient and useful for
FMM algorithm.
3.4 Numerical results
The algorithm outlined above has been implemented and extensively tested.
We have used LAPACK for solving the eigenvalue problem, where the reverse com-
munication function is ideally suitable for the inclusion of FMM codes for matrix-
vector products. Figure 3.12 compares the speed of FMM algorithm and the direct
summation scheme. It can been seen that the matching point is about N = 3000
where FMM starts to be faster than direct summation. Figure 3.13 illustrates the
error versus the truncation number p. Smaller p provides faster speed but compro-
mises the accuracy. So it is always a tradeofi in this situation. In our computations,
p = 5 is used in most of the cases. The number of levels or the number of the
sources in the flnest box is also a parameter which needs to be optimized to achieve
higher speed in FMM computation. Figure 3.14 tells us that the relation between
CPU time and group parameters s is nonlinear.
Table 3.1 presents the comparison between the results computed by the using
outlined technique and reported in [13] for the three silver spherical nanoparticle
system (see Figure 3.15). The experimentally measured dispersion relation for silver
published in [68] has been used in calculations. This table reveals a very good
agreement between our results and experimental results reported in [83].
The testing of the method has been performed for the more complicated sys-
76
tem of 20nm gold nanospheres shown in Figure 3.16. The separation between the
central sphere and surrounding spheres is 7nm. This structure is used to model the
aggregation of gold nanospheres due to DNA hybridization [84]. The low density
DNA solution is treated in computations as water. The calculations show that the
resonance wavelength for a single nanosphere is 516nm, while the resonance wave-
length for the group of 7 nanospheres shifts to 554nm. These results agree with
the experimental results reported in [84] where it is reported that for single gold
nanospheres the resonances occur at 525nm, while the resonances for the aggrega-
tion of gold nanospheres due to DNA hybridization occur at 560 nm. The above
computations have been performed on a 3.0GHz Dell-workstation.
Table 3.2 presents the comparison between FMM and conventional direct cal-
culations for the particle conflguration shown in Figure 3.16. Table II clearly reveals
the e?ciency of the fast multipole method, which is the method of choice for the
large size problems.
Table 3.1: Comparison of the resonance wavelengths for three silver spherical
nanoparticle system.
Case 1 Case 2
Conflguration ri;i+1=ri = 1=3, ri;i+1=ri = 1=3
ri;i+1=ri = 0:3, ri;i+1=ri = 0:6
Resonance wavelength [13] 369 nm 382 nm
Resonance wavelength (Computational) 372 nm 387 nm
77
102 103 104 105 10610
?2
10?1
100
101
102
103
Number of points N
CPU time (s)
Direct summation 
FMM algorithm 
y = bx 
y = cx2 
Figure3.12: TheCPUtimeinsecondsrequiredforcomputationofproblemsofsizeN
using direct method and FMM. FMM computations using truncation number p = 5.
Computation is performed on Dell workstation with Pentium 3.0 GHz processor, 4.0
GB RAM.
78
0 50 100 150 200 25010
?15
10?10
10?5
100
Truncation numbers p2
Absolute error
Figure 3.13: Absolute error vs truncation numbers p2 for the case N = 14496.
Computation is performed on Dell workstation with Pentium 3.0 GHz processor,
4.0 GB RAM.
79
10 20 30 40 50 60 70 80 90 10010
1
102
Number of points in smallest box 
CPU time (s)
Figure 3.14: CPU times vs group parameters s (number of points in the smallest
box) for N = 14496. Computation is performed on Dell workstation with Pentium
3.0 GHz processor, 4.0 GB RAM.
80
Figure 3.15: Three-sphere system for the nano-lens proposed in [83].



Figure 3.16: The schematic of a 7-sphere cluster.
81
Table 3.2: Computation times in seconds.
Method/mesh 8690 33600 134400
FMM 26 118 680
Direct summation 315 14819 237920*
*This computation time is estimated by extrapolation.
3.5 Multipole expansion method for spherical nanoparticles
In this section, a meshless method for the analysis of plasmon resonances in
multiple spherical nanoparticles is presented. It is shown that for spherical nanopar-
ticles, the eigenvalue problem deflned by (2.23) can be analytically simplifled. For
the case of coaxial spheres, the problem is appreciably simplifled due to its sym-
metry. Numerical results and computational e?ciency are presented and compared
with boundary element techniques.
3.5.1 The dual formulation for multiple spheres
+iV
0?? =? ( )??? ++=
n

iS
?V
i?1S
2S
3S
NS
Figure 3.17: Multiple spherical nanoparticles located in proximity.
82
Consider N spherical nanoparticles located in proximity to one another with
the same permittivity ?+(!) (Fig. 3.17). The source-free electric displacement D
satisfles everywhere
r?D = 0; (3.55)
r?D = 0; (3.56)
Then a scaler potential ' can be introduced
D = ?r'; (3.57)
such that ' is a solution to the Laplace equation
r2' = 0; (3.58)
inside (V + = V +1 [???[ V +N ) and outside (V ?) of the union of the nanoparticles
subject to the following interface boundary conditions on S, where S = S1[???[SN:
'+
?+(!) =
'?
?0 ; (3.59)
@'+
@n =
@'?
@n ; (3.60)
as well as the condition at inflnity
'(1) = 0: (3.61)
The potential of this sourceless fleld can be represented as a potential of double
layer of electric charges ?(x) distributed over S:
'(y) = 14??
0
NX
j=1
I
S
?(x)
jx?yjdSx; (3.62)
83
where y and x are points on surface S. In other words, a double layer of electric
charges ?(x) distributed on S may create the same electric displacement in the
free space as the resonant source-free electric displacement that may exist in the
presence of the nanoparticles. This will be the case if the potential (3.62) satisfles
all the conditions of boundary value problem (3.58){(3.61). It is apparent that the
potential (3.62) of the surface charges ?(x) satisfles the Laplace equation in V + and
V ?, decays to zero at inflnity and is continuous across S. To satisfy the boundary
condition (3.60), we recall that the normal derivatives of a single layer potential are
given by the formulas[70, 71]:
@'?
@n (y) = ?
?(y)
2?0 +
1
4??0
I
S
?(x)(x?y)?nxjx?yj3 dSx: (3.63)
where nx is the outward normal at point x. By substituting (3.9) into the boundary
condition (3.60), we arrive at the homogeneous boundary integral equation:
?(i)(y) = ?2?
NX
j=1
I
Sj
?(j)(x)(y?x)?nxjy?xj3 dSx; (3.64)
where
? = ?+(!)??0?
+(!)+?0
: (3.65)
3.5.2 Multipole expansion method
It is natural to apply the standard boundary element method to solve the
integral equation (3.64). For multiple spherical nanoparticles, however, a semi-
analytical method can be developed to simplify the problem which provides much
less complexity and computational costs in comparison with BEM. The main idea
84
of the approach can be stated as follows. Since the unknown functions ? in equation
(3.64) are deflned on the surfaces of the spheres, they can be represented by spherical
harmonics expansions. On the other hand, the kernel of the integral equation is
1=r-type, so it can be also expanded in terms of spherical harmonics via addition
theorem. Then by applying the orthogonality relations for spherical harmonics the
integrals can be evaluated analytically. If the \integration" point x and \evaluation"
point y are located on the same sphere, the spherical harmonics expansions possess
the same basis and the corresponding integrals can be evaluated directly. If the
\integration" point x and \evaluation" point y are located on difierent spheres,
the translation formulas for spherical harmonics have to be applied to move the
expansion center of \integration" sphere to the \evaluation" sphere such that they
have the same basis function. Finally, the linear equation for expansion coe?cients
can be obtained.
We assume that the total dimension of the system is appreciably less than
light wavelength (i.e., in quasi-static regime). The unknown surface function ?(z)
can be expanded into spherical harmonics
?(i)(z) =
1X
l=0
lX
m=?l
fi(i)lmYlm( dz?ci); (3.66)
where fi(i)lm are the expansion coe?cients, l = 0;1;2:::;m = ?l;:::;l, z is an
arbitrary point on the surface Si, ci is the center of Si, and dz?ci is a unit vector
deflned as z?cijz?cij.
If the \evaluation" point y and \integration" point x are located on the same
sphere Si which is the case for the ith term of the right hand side of integral equation
85









?

Figure 3.18: Points x and y are located at same sphere.
(3.64), the integral can be simplifled as follows. From Fig. 3.18, it is easy to see
that:
(x?y)?nx
jx?yj2 =
cosfl
jx?yj =
1
2ri; (3.67)
Using formulas (3.66) and (3.67), The ith term mentioned above can be written as:
I
Si
?(i)(x)(y?x)?nxjy?xj3 dSx = 12r
i
1X
l=0
lX
m=?l
fi(i)lm
I
Si
Ylm( dx?ci)
jy?xj dSx; (3.68)
By applying the addition theorem, the term 1jy?xj can be expanded into spherical
harmonics with center ci,
I
si
Ylm( dx?ci)
jy?xj dSx =
1X
l0=0
l0X
m0=?l0
4?
2l0 +1
I
si
Ylm( dx?ci)Y ?l0m0( dx?ci)Ylm(
dy?ci)
ri r
2
i sin d d`:
(3.69)
According to the orthogonality relations for spherical harmonics,
Z 2?
0
d`
Z ?
0
sin d Y ?l0m0( ;?)Ylm( ;?) = ?l0l?m0m; (3.70)
and obtain the flnal result for the ith term:
I
si
?(i)(x)(y?x)?nxjy?xj3 dSx = 2?
1X
l=0
lX
m=?l
fi(i)lmYlm(
dy?ci)
2l +1 ; (3.71)
86

 










Figure 3.19: Points x and y are located at difierent spheres.
If the \evaluation" point y and \integration" point x are located on the difier-
ent sphere Si and Sj, the integral can be simplifled as follows. Recall (r?)?n =@?@n
and use formula (3.66), The jth term of the right hand side of integral equation
(3.64) becomes
I
sj
?(j)(x)(y?x)?nxjy?xj3 dSx = ?
I
si
2
4
1X
l=0
lX
m=?l
fi(j)lmYlm( dx?cj)
3
5 @
@nx
? 1
jy?xj
!
dSx;
(3.72)
By using the addition theorem and @@nx = @@rj, we get
@
@nx
? 1
jy?xj
!
=
1X
l0=0
l0X
m0=?l0
4?
2l0 +1
@
@rj
0
@ jx?cjj
l0
jy?cjjl0+1
1
AY ?l0m0( dx?cj)Yl0m0( dy?cj)
=
1X
l0=0
l0X
m0=?l0
4?l0rl0?1j
2l0 +1
Yl0m0( dy?cj)
jy?cjjl0+1 Y
?
l0m0( dx?cj); (3.73)
Then according to the orthogonality relations for spherical harmonics, the integral
87
in (3.19) can be evaluated as
I
sj
?(j)(x)(y?x)?nxjy?xj3 dSx = ?4?
1X
l=0
lX
m=?l
fi(j)lm lr
l+1
j
2l +1
?Y
lm( dy?cj)
jy?cjjl+1
!
: (3.74)
In order to have the same basis functions on both sides of equation (3.64), we
translate the expansions from the center of \integration" sphere Sj to the center
of \observation" sphere Si. The re-expansion coe?cients  (j)l0m0 can be related to
the original expansion coe?cients fi(j)lm by using the \multipole-to-local" translation
operator and expressed in the following form
 jl0m0 =
1X
l=0
lX
m=?l
ijm0?mj?jm0j?jmjAml Am0l0
(?1)lAm?m0l0+l
Y m?m0l0+l (bt)
jtjl0+l+1
lrl+1j rli
2l +1fi
(j)
lm; (3.75)
where vector t is deflned in Fig. 3.19 and Aml is given by (3.16). To keep things
simpler, denote
 jl0m0 = K(j)l0m0;lmfi(j)lm; (3.76)
where K(j)l0m0;lm is the operator to perform the translation in equation (3.75). By
using this translation operator, equation (3.74) becomes:
I
Sj
?j(x)(y?x)?nxjy?xj3 dSx = ?4?
1X
l0=0
l0X
m0=?l0
K(j)l0m0;lmfi(j)lmYlm( dy?ci): (3.77)
By Substituting formulas (3.71) and (3.77) into (3.64) and equating the coe?cients
for same terms, we get
(1? ? 12l +1)fi(i)lm ?
NX
j=i;j6=i
fi(i)lmK(j)l0m0;lm = 0: (3.78)
We can write above equation for each sphere, then we get the simplifled eigenvalue
problem:
AX = flX; (3.79)
88
where
X =
h
fi(1)lm;???fi(i)lm;???fi(N)lm
iT
; (3.80)
fl = 1?; (3.81)
Aji = f
h 1
2l+1
i
; j = i;
K(j)l0m0;lm; j 6= i;
(3.82)
where i;j = 1;2;:::;N. It should be noted that Aji is a submatrix for A. When
j = i, it is a diagonal matrix, and when j 6= i, it is a dense matrix determined by
the translation operator.
In numerical simulations, the inflnite series of spherical harmonic has to be
truncated. If the truncation number is p, the total length of the series is p2. For
N-sphere system, the dimension for matrix A is Np2, and
?i(z) ?
p?1X
l=0
lX
m=?1
fi(i)lmYlm(z?ci): (3.83)
The truncation number p is not only a function of the error but also afiected by
the geometric conflguration. Smaller distances between spheres require a larger
truncation number to keep the accuracy. This phenomenon is shown in Fig.5. In
our simulation, we truncate as p ? pmax = 5rmax=dmin, where rmax is the radius of
the largest sphere in the system, dmin is the smallest gap between spheres. We have
checked that the numerical error associated with such a truncation is negligible.
If all the centers of spheres stay along a line, this is the coaxial case. Due to
the axial symmetry, the multipole to local translation operator can be decomposed
89
2 3 4 5 6 7 8 9 10 11 121.9
2
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
gap = 1 
gap = 0.5 
gap = 0.2 
Figure 3.20: Truncation number as a function of the gap between two spheres.
90
Table 3.3: Computed eigenvalue for two spheres with truncation number p = 10.
Mode Eigenvalue Mode Eigenvalue Mode Eigenvalue Mode Eigenvalue
1 2.6190 11 5.0276 21 6.9635 31 8.9070
2 2.8049 12 5.0742 22 7.0376 32 8.9330
3 2.8049 13 5.1823 23 7.1145 33 8.9330
4 3.2189 14 5.3208 24 7.1146 34 8.9845
5 3.2189 15 5.3209 25 7.1314 35 9.0198
6 3.4245 16 5.6506 26 7.1315 36 9.0199
7 4.7415 17 6.8485 27 7.1635 37 9.0391
8 4.7435 18 6.8837 28 7.3933 38 9.0594
9 4.7435 19 6.8837 29 7.3933 39 9.1431
10 4.9275 20 6.8931 30 7.6609 40 9.1433
and therefore the flnal matrix can be simplifled. The operator can be expressed as
following:
?
SjR
?
(t) = T(m;n;l;k)P
m?k
l+n (cosfi)e
imfl
?l+n+1 (3.84)
where m;n;l;k are indexes, and T is a known function of indexes. (?;fi;fl) represents
the translation vector t in spherical coordinate. In the coaxial case, because fl = 0,
the elements for difierent m are decoupled, and the matrix of can be decomposed
according to index m. For example, with truncation number p = 10, the size of the
translation matrix is 100. If we solve the problem for the two-sphere system, the
total size is 200. This matrix could be decomposed according to difierent m and
solved separately where the size of largest submatrix is 20. The results are shown
in Table 3.3 and Table 3.4.
91
Table 3.4: Computed eigenvalue for two spheres with truncation number p = 10
using submatrix.
m=0 2.6190 3.4245 4.7415 5.6506 7.1145 7.6609 ...
m=-1 2.8049 3.2189 4.7435 5.3208 7.0376 7.3933 ...
m=1 2.8049 3.2189 4.7435 5.3209 6.9635 7.3933 ...
m=-2 4.9275 5.0742 6.8931 7.1635 ...
m=2 5.0274 5.1823 6.8837 7.1315 ...
m=-3 6.8837 7.1314 ...
m=3 6.8485 7.1145 ...
... ... ... ... ... ... ... ...
3.5.3 Numerical results
The multipole expansion method has been numerically implemented and ex-
tensivelytested. Ithasprovedtobeaccurateande?cient. Theresultsareillustrated
by the examples presented below.
For single sphere, the results have been shown to be analytical and agree
with classical Mie theory. For two-sphere system, it always coaxial and therefore
can be further decomposed according to m. We computed the case with geometric
parameters: radius of spheres r1 = 1;r2 = 0:5, the distance between the surfaces
of the spheres d = 0:5. For the three-sphere system, we tested two difierent cases:
one is coaxial with parameters: r1 = 1;r2 = 0:5;r3 = 1;d12 = d21 = 0:5; the
other one, the three spheres are equal and form a equilateral triangle with r1 =
r2 = r3 = 1;d12 = d23 = d34 = 1. The results as well as the computation costs
are compared with BEM results and shown in Table 3.5 and Table 3.6. The results
of the two methods agree with each other very well and the multipole expansion
method is obviously more e?cient. In fact, this semi-analytical method also has
92
Table 3.5: Comparison of numerical results of BEM and MEM for two-nanosphere
and three-nanosphere conflgurations.
two-sphere three-sphere
Mode BEM MEM Mode BEM MEM
1 2.4959 2.5161 1 2.2659 2.3453
2 2.7378 2.7290 2 2.6681 2.6594
3 2.7379 2.7290 3 2.6681 2.6594
4 3.2600 3.3391 4 2.9028 2.9427
5 3.2603 3.3391 5 2.9028 2.9427
6 3.4666 3.5652 6 3.1300 3.1262
Table 3.6: Comparison of computation costs of BEM and MEM for two-nanosphere
and three-nanosphere conflgurations.
Method Number of sphere matrix size CPU time(s)
MEM 2 200 0.1
BEM 2 2560 3.8 (with ARPACK)
MEM 3 300 0.25
BEM 3 3840 8.7 (with ARPACK)
the advantage for assembly of the matrix, because it only requires to build up the
translation matrix. Each non-diagonal submatrix has the same structure; it repeats
only with difierent translation vectors.
93
Chapter 4
Extinction Cross Sections of Nanoparticles
In this chapter, the boundary integral equation technique (developed in chapter 2)
is extended to the calculation of extinction cross sections (ECS) of nanoparticles.
ECS is an important measure of optical properties of nanoparticles which reveals
the strength and the sharpness of the plasmon resoance modes. The concept of
ECS and the optical theorem are brie y reviewed in section 4.1. The computation
algorithm of ECS as well as the numerical results are presented in sections 4.2 and
4.3. It is shown that the scattered fleld can be derived through the solution of an
inhomogeneous integral equation, and the extinction cross section is then obtained
by employing the optical theorem. Finally, the outline technique is modifled to
take into account the phase shift between adjacent nanoparticles for the analysis of
nanoparticle-structured plasmonic waveguides of light in section 4.4.
94
4.1 What is the extinction cross section ?
4.1.1 Concept of extinction cross section
The extinction cross section (ECS) of a nanoparticle represents the total power
loss from the incident wave due to the scattering and absorption by the nanoparticle.
The ECS curve clearly reveals the plasmon resonances by peaks around resonance
frequencies and itprovidesimportantinformation suchasthe strength and sharpness
of the resonance. The mathematical deflnition of ECS can be brie y demonstrated
as follows [85, 86].
Figure 4.1: Scattering diagram of a nanoparticle.
Consider the scattering geometry in Figure 4.1, a plane wave incident on a
scatterer with arbitrary shape. The fleld (E?) at any point in the lossless medium
with permittivity ? and permeability ? surrounding the scatterer may be represented
95
as the sum of the incident fleld (Ei) and scattered fleld (Es).
E? = Ei +Es; (4.1)
H? = Hi +Hs: (4.2)
Let the permittivity of the scatterer (nanoparticle) be complex valued with real and
imaginary parts
?+(!) = ?0+(!)+i?00+(!): (4.3)
The electromagnetic flelds inside the particle (V +) obey the equations:
r?E+ = j!?H+; (4.4)
r?H+ = ?j!?+E+; (4.5)
and are subject to the following boundary conditions:
n?E+ = n?(Ei +Es); (4.6)
n?H+ = n?(Hi +Hs): (4.7)
The time-averaged power absorbed by the nanoparticle is
Wa = 12
Z
V +
!?00+
flfl
flE+
flfl
fl2 dV: (4.8)
It is easy to show that:
r?
?
E+ ?H+?
?
= H+??r?E+?E+?r?H+? = j!?
flfl
flH+
flfl
fl2?j!??+
flfl
flE+
flfl
fl2 ; (4.9)
and taking the real part of (4.9) gives:
r?
h
Re
?
E+ ?H+?
?i
= ?!?00+
flfl
flE+
flfl
fl2 : (4.10)
96
Therefore the absorbed power by the nanoparticle can be expressed as:
Wa = ?12
Z
V +
r?
h
Re
?
E+ ?H+?
?i
dV = ?12
Z
S
n?Re
?
E+ ?H+?
?
dS: (4.11)
Substitute boundary conditions (4.6) and (4.7) into formula (4.11), it is easy to see
that:
Wa = ?12
Z
S
n?Re[(Ei +Es)?(H?i +H?s)]dS = ?
Z
S
n?(Si +S0 +Ss)dS; (4.12)
where Si;S0;Ss are deflned as:
Si = 12Re(Ei ?H?i); (4.13)
S0 = 12Re[(Ei ?H?s)+(Es ?H?i)]; (4.14)
Ss = 12Re(Es ?H?s): (4.15)
Since r? Si = 0 is true everywhere in V + and V ?, the incident power Si should
satisfy
Z
S
n?SidS = 0; (4.16)
and because there is no source between boundary S and inflnity, the scattered power
can be evaluated as an integral over S:
Ws =
Z
S
n?SsdS: (4.17)
Then the absorbed power now can be written in the following form:
Wa = ?
Z
S
n?S0dS ?Ws; (4.18)
and the total power loss due to the scatterer (scattering and absorption) is:
Wext = Wa +Ws = ?
Z
S
n?S0dS: (4.19)
97
Thus the extinction cross section (ECS) is deflned as the ratio of total power Wext
to the incident power per unit area:
 ext = WextI
i
; (4.20)
where Ii is the incident radiation. Similarly, we can deflne the scattering cross
section (SCS) and the absorption cross section (ACS):
 sca = WsI
i
; (4.21)
 abs = WaI
i
; (4.22)
and it is apparent:
 ext =  sca + abs: (4.23)
It is instrumental to point out that for small metallic particles, the losses of incident
radiation is dominated by absorption, so the above express can be written as:
 ext ?  abs: (4.24)
4.1.2 The optical theorem
The above deflnition can be used to compute the ECS of a nanoparticle. An-
other convenient way, however, is to employ the so-called optical theorem. The
optical theorem expresses a very curious fact: the extinction depends only on the
scattering amplitude in the forward direction. The proof of this theorem can be
found in [39, 87] and will be omitted here. In the following, only those key formulas
are summarized for the calculation of ECS by using the optical theorem.
98
Let the incident wave be linearly polarized with polarization vector bei:
Ei = beiE0eiki?r; (4.25)
Hi =
?b
ki ?bei
? E0
? e
iki?r; (4.26)
Bi = 1ckbki ?Ei: (4.27)
The far fleld can be related to the scattering amplitude F
?b
ks;bki
?
as:
F
?b
ks;bki
?
= limr!1
?eikr
r
!?1
Es (r); (4.28)
where ki is the incident wave vector, ks is the scattering wave vector in the direction
of observation, Es is the scattered fleld in far-fleld zone. Then the optical theorem
tells us:
 ext = 4?k Im
h
be?i ?f
?b
ki;bki
?i
; (4.29)
where k is the wave number, f
?b
ki;bki
?
is the normalized forward direction scattering
amplitude and the normalized scattering amplitude is deflned as:
f
?b
ks;bki
?
= F
?b
ks;bki
?
E0 : (4.30)
According to (4.28){(4.30), the entire problem of the computation of ECS is reduced
to the evaluation of the forward scattering amplitude F, which is related to the
scattered far fleld Es.
4.1.3 Measurement of ECS
In practice, the ECS can be directly measured (see Fig.2) in terms of the
following expression [39]:
U = Ii[A(D)? ext]; (4.31)
99


	



Figure 4.2: Schematic illustration for the measurement of extinction cross sections.
where Ii is the incident radiation, A(D) is the area of the detector and U is the
power incident on the detector. We measure U with and without the particle inter-
posed between source and detector. Because  ext is inherently positive, the efiect
of the particle is to reduce the detector area by  ext. This method of measurement
illustrates the interpretation of  ext as an area.
4.2 Computation of ECS
4.2.1 Near fleld computation
Consider an incident plane wave scattered by nanoparticle of arbitrary shape
with permittivity ?+(!) (Fig. 4.1). We shall start with Maxwell?s equations written
in the mathematical form that explicitly re ects the smallness of particle dimensions
in comparison with the free-space wavelength. For this reason, we introduce the
100
scaled vectors of scattered and incident flelds
e = ?120E; h = ?120H; ein = ?120Ein (4.32)
as well as spatial coordinate scaled by the diameter d of the object. This leads to
the following boundary value problem for scattered flelds e? and h?:
r?e+ = ?jflh+; r?h+ = jfl?+?
0
e+ +jfl
 ?
+
?0 ?1
?
e+in; (4.33)
r?e+ = 0; r?h+ = 0; (4.34)
r?e? = ?jflh?; r?h? = jfle?; (4.35)
r?e? = 0; r?h?;= 0 (4.36)
n?
?
e+ ?e?
?
= 0; n?
?
h+ ?h?
?
= 0; (4.37)
n?
?
?+e+ ??0e?
?
= (?0 ??+)n?e+in; n?
?
h+ ?h?
?
= 0; (4.38)
where superscripts \+" and \-" are used for physical quantities inside (V +) and
outside (V ?) the dielectric object, respectively, n is a outward unit normal to S,
while
fl = !p?0?0d: (4.39)
In the case when the free-space wavelength is large in comparison with nanoparticle
dimension, fl can be treated as a small parameter and solution of the boundary
value problem (4.33)-(4.38) can be expanded in terms of fl
e? = e?0 +fle?1 +fl2e?2 +???; (4.40)
h? = h?0 +flh?1 +fl2h?2 +???: (4.41)
101
An incident fleld is usually treated as independent of the characteristic dimension of
the system under consideration. For this reason, it can be regarded as independent
of fl. By substituting formulas (4.40){(4.41) into equations (4.33){(4.36) as well
as boundary conditions (4.37){(4.38) and equating terms of equal powers of fl, we
obtain the boundary value problems for e?k and h?k (k = 0;1;2:::).
For zero-order terms, these boundary value problems can be written in terms
of E?0 = ??120 e?0 and H?0 = ??120 h?0 as follows:
r?E?0 = 0; r?E?0 = 0; (4.42)
n?
?
E+0 ?E?0
?
= 0; n?
?
?+E+0 ??0E?0
?
= (?0 ??+)n?E+in; (4.43)
and
r?H?0 = 0; r?H?0 = 0; (4.44)
n?
?
H+0 ?H?0
?
= 0; n?
?
H+0 ?H?0
?
= 0: (4.45)
From (4.44)-(4.45), it is apparent that
h?0 = 0: (4.46)
The electric potential ? can be introduced for the electric fleld E0 and this potential
can be represented as an electric potential of a single layer of electric charge  
distributed over the boundary S of the particle. It is clear that the electric fleld
of surface charges is curl and divergence free in V + and V ? and satisfles the flrst
boundary condition in (4.43). Next we recall the properties of single layer potential
(2.20). By substituting (2.20) into the second boundary condition in (4.43), after
102
simple transformation we arrive at the following inhomogeneous boundary integral
equation:
 0 (Q)? ?2?
I
S
 0 (M) nQ ?rMQr3
MQ
dSM = 2?0?nQ?e+in (Q); (4.47)
where
? = ?+(!)??0?
+(!)+?0
: (4.48)
Thus, the scattered electric flelds can be computed (in zero-order) as follows:
e?0 = ?r?0 = ?r
" 1
4??0
I
S
 0 (M)
rMQ dSM
#
: (4.49)
The inhomogeneous boundary integral equation (4.47) for the calculations of electro-
static and scattering problem with negative ? was extensively studied in publications
of Mayergoyz [54, 55] and introduced for the analysis of plasmon resonance in [20].
From equations (4.33){(4.38), (4.40){(4.41) and (4.46), the boundary value problem
for the flrst-order corrections can be derived and appears as the following:
r?e?1 = 0; r?e?1 = 0; (4.50)
n?
?
e+1 ?e?1
?
= 0; n?
?
?+e+1 ??0e?1
?
= 0; (4.51)
r?h+1 = j?
+
?0 e
+
0 +j
 ?
+
?0 ?1
?
e+in; r?h?1 = je?0 ; r?h?1 = 0; (4.52)
n?
?
h+1 ?h?1
?
= 0; n?
?
h+1 ?h?1
?
= 0: (4.53)
From equations (4.50){(4.51), it is clear that
e?1 = 0: (4.54)
Next, we proceed to the solution of boundary value problem (4.52)-(4.53). Terms
j?+?0 e+0 +j
??
+
?0 ?1
?
e+in and je?0 in the flrst two equation of (4.52) can be interpreted
103
as current sources and the solution of boundary value problem (4.52)-(4.53) can be
written in the integral form
h?1 (Q) = j
?+
?0
4?
Z
V +
rMQ
r3MQ ?e
+
0 (M)dVM +
j
4?
Z
V?
rMQ
r3MQ ?e
?
0 (M)dVM
+ j
??
+
?0 ?1
?
4?
Z
V +
rMQ
r3MQ ?e
+
in (M)dVM: (4.55)
The last expression can be appreciably simplifled and reduced to an integral over
boundary S. Indeed, by using the fact that r ? e?0 = 0 and by employing the
following formulas for vector analysis:
r?(gF) = (rg)?F+g(r?F); (4.56)
Z
V
dV (r?A) =
I
S
dS?A; (4.57)
after simple transformations we arrive at
h?1 (Q) = ?j
??
+
?0 ?1
?
4?
8<
:
I
S
nM ?
h
e+0 (M)+e+in (M)
i
rMQ dSM
9=
;: (4.58)
Now we proceed to the discussion of second order terms in expansions of e+ and h+.
From equations (4.33){(4.38), (4.40){(4.41) and (4.54){(4.55), the boundary value
problem for the second-order terms e+2 and h+2 can be derived, respectively:
r?e?2 = ?jh?1 ; r?e?2 = 0; (4.59)
n?
?
?+e+2 ??0e?2
?
= 0; n?
?
e+2 ?e?2
?
= 0; (4.60)
and
r?h?2 = 0; r?h?2 = 0; (4.61)
n?
?
h+2 ?h?2
?
= 0; n?
?
h+2 ?h?2
?
= 0: (4.62)
104
From formulas (4.61){(4.62), it is apparent that
h?2 = 0: (4.63)
In order to solve the boundary value problem (4.59){(4.60), we shall split the electric
fleld e?2 into two distinct components
e?2 = e?2 + ee?2 ; (4.64)
that satisfy the following boundary value problems, respectively:
r?e?2 = ?jh?1 ; r?e?2 = 0; (4.65)
n?
?
e+2 ?e?2
?
= 0; n?
?
e+2 ?fe ?2
?
= 0; (4.66)
and
r?ee?2 = 0; r?ee?2 = 0; (4.67)
n?
?
?+ee+2 ??0ee?2
?
= (?+ ??0)n?e+2 ; n?
?e
e+2 ?ee?2
?
= 0: (4.68)
By using the same line of reasoning as in section 2.2, it can be shown that the
solution of the boundary value problem (4.65)-(4.66) can be written as follows
e?2 = ?
??
+
?0 ?1
?
8?
I
S
nM ?
h
e+0 (M)+e+in (M)
i
?rMQ
rMQ dSM: (4.69)
Next, we proceed to the solution of the boundary value problem (4.67){(4.68).
Again, the single layer charge distribution can be introduced to construct the elec-
tric fleld ee?2 . Then by using the same line of reasoning in the derivation of equation
(2.51), we obtain the following inhomogeneous integral equation for  2(M):
 2 (Q) ? ?2?
I
S
 2 (M) nQ ?rMQr3
MQ
dSM
= 2?0?+?
??
+
?0 ?1
?
8? ?
I
S
nQ ?
h
e+0 (M)+e+in (M)
i
?rMQ
rMQ dSM: (4.70)
105
Thus, the second-order term for scattered electric fleld e+2 can be written as:
e?2 = ?r
" 1
4??0
I  
2 (M)
rMQ dSM
#
?
??+
?0 ?1
?
8?
I
S
nM ?
h
e+0 (M)+e+in (M)
i
?rMQ
rMQ dSM:
(4.71)
Finally, the scattered near fleld is given by the formula
E+ ? ??120
?
e+0 +fl2e+2
?
: (4.72)
4.2.2 Far fleld computation
As discussed in section 4.1, the key step to obtain ECS is the computation of
scattered far flelds. In the following, we demonstrate the technique to compute this
quantity according to the obtained scattered near fleld. The far fleld is the solution
of the following boundary value problem:
r?E? = ?j!?0H?; (4.73)
r?H+ = j!?0E+ +j!(???0)
?
E+ +E+in
?
; (4.74)
r?H? = j!?0E?; (4.75)
n?
?
E+ ?E?
?
= 0; n?
?
H+ ?H?
?
= 0; (4.76)
where E? = ??120 e?, H? = ??120 h?. In this case, the near fleld can be interpreted
as the induced polarization current in V +; in other words, the flctitious polariza-
tion current can be determined by the known near fleld. Then, the far fleld can
be computed as the flelds created by the polarization current in free-space. The
polarization current is deflned as
iP = @P@t ; (4.77)
106
where P is the polarization
P = (?+ ??0)E+total: (4.78)
For time-harmonic flelds, it is easy to see:
iP = j!(?+ ??0)
?
E+ +E+in
?
: (4.79)
On the other hand, the scattered far fleld can be expressed in terms of localized
polarization current as
h(Q) = r?A(Q); (4.80)
and
A(Q) = ?04?
?Z
V +
iP (M)dVM
??e?ikrMQ
rMQ
!
: (4.81)
By using formula (4.56), from (4.80)-(4.81) we get
h?(Q) = 14?
Z
V +
iP (M)?rQ
?e?ikr
MQ
rMQ
!
dVM: (4.82)
Note that
rQ
?e?ikr
r
!
= ?rMQ
 
ik + 1r
? e?ikr
r ; (4.83)
and for the scattered fleld in the far fleld zone, we can drop the 1=r term; therefore,
equation (4.82) is reduced to
H?(Q) ? ik4?
Z
V +
[brMQ ?iP (M)] e
?ikrMQ
rMQ dVM: (4.84)
The last expression can be further simplifled since the term brMQe?ikrMQr
MQ
does not
depend on point M explicitly when the observation is made from the far fleld
H?s (Q) ? ck
2
4?
e?ikrMQ
rMQ brMQ ?
?Z
V +
P(M)dVM
?
: (4.85)
107
Since the electromagnetic flelds in the far fleld zone can be always considered as
plane waves, the electric fleld and magnetic fleld has the following simple relation
E?s (Q) = Z0H?s (Q)?k; (4.86)
where Z0 is the intrinsic impedance of free-space. So we get the flnal expression for
electrical far fleld:
E?s (Q) = ck
2Z0
4?
e?ikrMQ
rMQ brMQ ?
?Z
V +
P(M)dVM
?
?brMQ: (4.87)
If we assume the polarization is along the x-axis, and propagation is along the z-axis,
according to (4.28), the forward scattering amplitude is
F(z;bz) = ck
2Z0
4? bz?
?Z
V +
P(M)dVM
?
?bz: (4.88)
Finally, by taking (4.88) into (4.29) we arrive following expression for the ECS:
 ext = ckZ0Im
"
bei ? (bz?
R
V + P(M)dVM)?bz
E0
#
: (4.89)
4.3 Numerical techniques to compute ECS
4.3.1 Discretization of the inhomogeneous integral equation
Let us partition S into N small pieces4Sj and rewrite integral equation (4.47)
as follows:
 (Qi)? ?2?
NeX
j=1
Z
?Sj
 (Mj) rMjQi ?nQir3
MjQi
dSj = 2?0?nQi?e+in: (4.90)
Now we integrate (4.90) over ?Si and exchange the order of the integral
Z
?Si
 (Qi)dSi? ?2?
NeX
j=1
Z
?Sj
 (Mj)
"Z
?Si
rMjQi ?nQi
r3MjQi dSi
#
dSj = 2?+?0?
Z
?Si
nQi?e+indSi:
(4.91)
108
By introducing the notation
!ij = f
rMjQi?nQi
r3MjQi ?Si if i 6= j;
2? ?PNej=1j6=i rMjQi?nQir3
MjQi
?Si if i = j;
(4.92)
the last formula can be presented as follows:
NeX
j=1
!ij
Z
?Sj
 0 (Mj)dSj ? 2??
Z
?Si
 (Qi)dSi = ?4??0
Z
?Si
nQi?e+in (Qi)dSi: (4.93)
Deflne:
Xi =
Z
?Si
 0 (Qi)dSi; (4.94)
Ii =
Z
?Si
nQi?e+in (Qi)dSi; (4.95)
we get the following linear system for surface charges:
NeX
j=1
!ijXj ? 2?? Xi = ?4??0Ii: (4.96)
4.3.2 Far fleld and dipole moment
Once the surface charge distribution is known, the far fleld can be computed
by using the algorithm introduced in previous section. However, there is a more
e?cient way for which the computation of far fleld can be reduced to the calculation
of dipole moment. In the following, we demonstrate how this can be done. According
to equation (4.85), the central problem is to compute the following integral:
IP =
Z
V +
P(M)dVM = (?+ ??0)
Z
V +
?
E+sc +E+in
?
dVM: (4.97)
Since the incident fleld is practically constant within the volume (because the free-
space wavelength is much larger then object dimension and we observe from far fleld
109
point), we have:
IPin = (?+ ??0)E+in(M0)V: (4.98)
This reduces the problem to how to calculate the polarization current due to scat-
tered flelds:
IPsc = (?+ ??0)??120
Z
V +
e+0 dVM: (4.99)
To this end, let us introduce the vector d and its scalar components as following:
d =
I
 (M)rMQdSM; (4.100)
dx =
I
 (M)xMdSM; dy =
I
 (M)yMdSM; dz =
I
 (M)zMdSM: (4.101)
On the other hand, according to formula (2.20), we have
 = ?0n?
?
E? ?E+
?
: (4.102)
By taking (4.102) into the flrst boundary condition in (4.38), after simple transfor-
mation, we arrive at
 = (?+ ??0)(n?E+ +n?Ein): (4.103)
Then the x-component of vector d can be computed as:
dx = (?+ ??0)
I ?
n?E+ +n?Ein
?
xMdSM = (?+ ??0)
Z
V
?
E+x +Einx
?
dVM:
(4.104)
Therefore, we have
d = fi
Z
V
?
E+ +Ein
?
dVM = IP: (4.105)
Then for the plane wave deflned in (4.24)-(4.26), according to (4.88) the ECS can
be represented in terms of d.
 ext = ckZ0Im
"
bx? (bz?d)?bzE
0
#
: (4.106)
110
This expression greatly simplifles the computation of ECS since the dipole moment
d is very easy to compute once the surface charge is known.
4.3.3 Numerical results
The numerical technique discussed above has been software implemented and
tested. The algorithm described above has been flrst tested for spherical particles
where exact analytical solutions are available (Mie theory).
 ext = ?a24xIm
? ?
p ??m
?p +2?m
!
; (4.107)
x = 2?? aN; (4.108)
where a is the radius of the sphere and N is the re ective index for surrounding
medium. We have computed the extinction cross section of nanospheres for difierent
materials. Fig. 4.3 and Fig. 4.4. are the results for Au and Ag, respectively. It is
apparent that the numerical results are very accurate. The computational results
of ECS for Au ellipsoidal nanoparticles with difierent aspect ratios are presented in
Figure 4.5. It can be seen that adjusting the aspect ratio of the ellipsoidal nanopar-
ticle is an efiective way to tune the plasmon resonances. Figure 4.6 demonstrates
the computational results for gold nano-rings placed on a dielectric substrate, illus-
trating that we have reasonable agreement with the experiment results. Finally, the
computation for semiconductor materials is also performed. The results for InSb
sphere are shown in Fig. 4.7.
111
400 450 500 550 600 650 7000
20
40
60
80
100
120
140
Wavelength (nm)
Extinction cross section (nm
2 )
computational results
Mie theory predictions
Figure 4.3: Computational results compared with the Mie theory predictions for the
extinction cross section of a single Au nanoparticle (diameter = 10 nm).
112
320 340 360 380 400 4200
500
1000
1500
2000
2500
Wavelength (nm)
Extinction cross?section (nm
2 )
computational results
Mie theory predictions
Figure 4.4: Computational results compared with the Mie theory predictions for the
extinction cross section of a single Ag nanosphere (diameter = 20 nm).
113
320 340 360 380 400 420 4400
500
1000
1500
2000
2500
3000
3500
Wavelength (nm)
Extinction cross section (nm
2 )
1 
2 
3 
4 
5 
Figure 4.5: Computational results of the extinction cross section of a single Ag
nanoellipsoid, with difierent aspect ratios: 1:1:0.8, 1:1:1. 1:1:1.2. 1:1:1.4. 1:1:1.6 for
curve 1, 2, 3, 4, 5, respectively.
114
800 1000 1200 1400 1600 18000
0.5
1
1.5
2
2.5
3x 10
5
Wavelength (nm)
Extinction cross?section (nm
2 )
14nm 
10nm 
9nm 
Figure 4.6: Computed extinction cross sections of Au nanorings placed on substrate,
the dimensions of these rings are provided in Table 2.3.
115
50 100 150 200 250 300 350 4000
100
200
300
400
500
600
700
Wavelength (nm)
extinction cross section (nm
2 )
Figure 4.7: Computational extinction cross section of a single InSb nanosphere
(diameter = 20 nm) placed on a glass substrate ? = 2:25?0.
116
4.4 Numerical analysis of plasmon waveguides of light
Plasmon waveguides of light have been the focus of considerable research lately
[10, 11, 12]. These plasmon waveguides consist of an array of metallic nanoparti-
cles with their resonance frequencies in the region of optical waveguiding. These
nanoparticle-structured waveguides hold the unique promise for light guiding and
bending at the nanoscale. For this reason, they are considered very instrumen-
tal in the emerging fleld of nanophotonics where the information transmission and
processing occur entirely at the optical level. In this section, the technique for the
calculation of resonance frequencies of nanoparticle-structured waveguides is pre-
sented. This technique is based on the boundary integral equation method and, in
this sense, it can be considered as the further extension of the technique for the
analysis of plasmon resonances in metallic nanoparticles developed in previous part
of this dissertation.
4.4.1 The numerical technique
To start the discussion, consider a chain of identical metallic nanoparticles
subject to incident optical radiation on its left edge (see Figure 4.8). At speciflc
frequencies, this optical incident radiation will excite plasmon resonances in the
flrst nanoparticle, which then through near-fleld coupling will induce plasmon reso-
nances in all nanoparticles that form the chain. This is, in a nutshell, the physical
mechanism of light plasmon waveguiding.
It is well known that plasmon resonances occur when the free-space wavelength
of light is large in comparison with particle dimensions. For this reason, the time-
117


	




	

Figure 4.8: Schematic of plasmonic waveguide of light.
harmonic electromagnetic flelds within each nanoparticle and around it vary almost
with the same phase. In other words, at any instant of time, the flelds locally
(around each particle) look like electrostatic flelds. By using this fact, the following
integral equation can be derived for each particle in the chain:
 p(Q)? ?2?
I
Sp
 p(M)rMQ ?nQr3
MQ
dSM = ?2??0nQ ?Eap; (4.109)
? = ?+(!)??0?
+(!)+?0
: (4.110)
Here, p is the nanoparticles number in the chain,  is the flctitious single layer of
charges distributed over boundary S, Eap is the electric fleld that has two distinct
parts (the incident fleld as well as the fleld scattered by all other nanoparticles), and
?+(!) and ?0 are the dielectric permittivities of nanoparticles and the surrounding
medium, respectively.
The dimensions of the waveguiding chain are usually comparable to the wave-
length of the incident radiation. As a result, the phase shift in time-harmonic
electromagnetic flelds scattered by difierent particles cannot be neglected. This
118
situation can be accounted for in equation (4.109) as follows:
 p(Q)? ?2?
I
Sp
 p(M)rMQ ?nQr3
MQ
dSM = ?2??0nQ ?(Ein +Escat); (4.111)
where Ein is the electric fleld of incident radiation, while Escat is the fleld scattered
by all other nanoparticles in the chain that can be computed by using the formula
Escat(Q) = 14??
0
I
Sp
 (M)
?
ik + 1r
MQ
! r
MQeikrMQ
r2MQ dSM;
~Sp = X
k6=p
Sk: (4.112)
By solving the coupled integral equations, the extinction cross-section and its de-
pendence on the frequency can be computed [8]. The frequency at which the ex-
tinction cross-section achieves its maximum value can be identifled as the resonance
frequency for which the light guiding is the most e?cient.
Next, webrie ydescribethediscretizationtechniqueforthesolutionofintegral
equation (4.111). In a similar procedure to what we did before, let us partition Sp
into N small pieces ?Sj, although the kernel in this case will be hybrid. Nevertheless,
equation (4.111) can be transformed to the following matrix form:
NX
j=1
!ijXj ? 2?? Xi +
MX
k=1;k6=p
NX
j=1
?km = ?4??0Ii; (4.113)
where
!i(M) =
Z
?Sj
rMQ ?nQ
r3MQ dSQ; (4.114)
Xi =
Z
?Sj
 (Q)dSQ; (4.115)
Ii =
Z
?Sj
nQ ?E(p)in nQ; (4.116)
?i =
Z
?Sj
nQ ?rMQ
?
ik + 1r
MQ
! r
MQeikrMQ
rMQ dSQ: (4.117)
119
The linear system deflned by (4.113) can be e?ciently solved by using the GM-
RES algorithm. The numerical technique based on discretization (4.113) has been
software implemented and extensively tested. This technique is illustrated by the
examples presented in the following section.
4.4.2 Numerical results
First, we present the results of computations and their comparison with avail-
able experimental data [11] for a chain of seven spherical gold nanoparticles of 50nm
diameter and separated by 25nm from one another. The comparison between the
resonance waveguide frequencies found experimentally in [11] and our computational
results for these frequencies is presented in the Table 4.1. This table reveals fairly
good agreement between computational and experimental results.
The extinction cross sections for longitudinal and transverse polarizations of
incident radiation are shown in Figure 4.9. It is apparent from this flgure that
the longitudinal mode has better transmission properties than the transverse mode;
this is not surprising because the longitudinal mode ofiers better conflnement of
plasmon resonant flelds. The silver and gold dispersion relations published in [9]
have been used in our calculations. Furthermore, we have performed computations
for the 11-gold-sphere waveguides with the difierent nanoparticles arrangements.
In these simulations, the nanoparticles of 50nm in diameter and 75nm center-to-
center distance have been used. For the chain structure (Fig. 4.10), the extinction
cross sections for the longitudinal mode are shown in Figure 4.11). By comparing
this flgure with Table 4.1, it can be observed that the resonance frequency of the
120
11-sphere chain is about 605nm, which is further shifted to the longer wavelength.
Figures 4.13 and 4.15 present the extinction cross-sections for the T-shape
(Fig. 4.12) and L-shape (Fig. 4.14) nanoparticle waveguides, respectively. For
the T-shape waveguide, the transverse polarization excitation is at the short end.
When the plasmon resonant flelds are coupled to the longer arm, the mode becomes
longitudinal. For the L-shape waveguide, the longitudinal polarization excitation is
at the end of the long arm. In these structures, the right angles between particle
arrays are responsible for the polarization transformations.
In this section, the novel numerical approach to the analysis of nanoparticle-
structured plasmon waveguides of light has been presented. The extinction cross-
sections and resonance (propagation) frequencies are computed for various geome-
tries of these waveguides and compared with available experimental data. The
computational results for resonance frequencies corresponding to difierent light po-
larizations are illustrated as well.
Table 4.1: Resonance frequencies for gold 7-sphere-chain waveguide.
Experimental result [11] Computational results
Transverse polarization 585 nm 534 nm
Longitudinal polarization 602 nm 572 nm
121
500 550 600 650 7000
500
1000
1500
2000
2500
3000
3500
4000
4500
Wavelength (nm)
Extunction cross?section
Transverse polarization
Longitudinal polarization
Figure 4.9: Computed extinction cross sections of gold 7-sphere-chain waveguide for
the transverse and longitudinal polarization modes (diameter = 50nm, center-to-
center distance = 25nm).
122
Figure 4.10: Schematic of gold 11-sphere-chain waveguide of light
500 550 600 650 700 7500
500
1000
1500
2000
2500
3000
3500
4000
4500
Wavelength (nm)
Extinction cross?sections
Figure 4.11: Computed extinction cross sections of gold 11-sphere-chain waveguide
for the longitudinal polarization mode (diameter = 50nm, center-to-center distance
= 25nm).
123
Figure 4.12: Schematic of gold L-shape 11-sphere waveguide of light
400 450 500 550 600 650 700400
600
800
1000
1200
1400
1600
1800
2000
2200
Wavelength (nm)
Extinction cross?sections
Figure 4.13: Computed extinction cross sections of gold L-shape 11-sphere
waveguide for the longitudinal polarization mode (diameter = 50nm, center-to-
center distance = 25nm).
124
Figure 4.14: Schematic of gold T-shape 11-sphere waveguide of light
480 500 520 540 560 580 600 620 640 660 6801000
2000
3000
4000
5000
6000
7000
8000
9000
Wavelength (nm)
Extinction cross?sections
Figure 4.15: Computed extinction cross sections of gold T-shape 11-sphere
waveguide for the longitudinal polarization mode (diameter = 50nm, center-to-
center distance = 25nm).
125
Chapter 5
Conclusion
Plasmon resonance in nanoparticles is a unique nanoscale phenomenon which
has numerous scientiflc and technological applications in such areas as near-fleld
microscopy, nano-lithography, surface enhanced Raman scattering, nanophotonics,
biosensors, optical data storage, etc. For this reason, the understanding of the
optical properties of nanoparticles holds both fundamental and practical signifl-
cance. While considerable experimental and theoretical work has been done in this
fleld, the accurate and e?cient modeling of plasmon resonances in three-dimensional
nanoparticles for arbitrary shape remains a challenge.
In this dissertation, we developed a robust and e?cient technique to fully char-
acterize the plasmon resonances and performed a comprehensive numerical analysis
of plasmon resonances in three-dimensional nanostructures. This surface integral
equation technique is based on the observations that plasmon resonances in nanopar-
ticles occur at speciflc frequencies for which the particle permittivity is negative and
the free-space wavelength of the radiation is large in comparison with particle di-
mensions. By using perturbation technique to Maxwell?s equations, the problem of
determining resonance frequencies is framed as a classical eigenvalue problem. It
is shown that the resonance values of dielectric permittivities as well as resonance
frequencies can be directly found through the solution of this eigenvalue problem
126
for three-dimensional nanoparticles.
Thisintegralequationmethodleadstofullypopulateddiscretizedmatrixequa-
tions that are computationally expensive to solve, especially when a large number
of particles are involved in the nanostructures. Since the fully populated matrices
are generated by integrals with 1=r-type kernel, this computational problem can
be appreciably alleviated by using the fast multipole method. This method greatly
speeds up the matrix-vector multiplications. By implementing the FMM algorithm
in the solution of the eigenvalue problem, the integral equation technique becomes
computationally very e?cient compared to other technique such as flnite-difierence
time-domain (FDTD) method.
We have extended this technique to the computation of extinction cross sec-
tions from which important information such as the strength and FWHM of plasmon
resonances can be obtained. The applications of plasmon resonances have also been
considered. We performed extensive numerical studies for metallic nanoshells and
plasmon waveguides of light, which are very promising for bio-sensing applications
and light guiding and bending under the difiraction limit, respectively. It can be seen
that the technique developed throughout this dissertation can be very instrumental
for the design of plasmon resonant nanoparticles and tailor their optical properties
for various applications.
To further advance the understanding of optical properties of nanoparticles
and pursue its applications in new promising technology, the following issues can be
foreseen in the near future:
(1) Time dynamic analysis of plasmon resonances. The temporal analysis of
127
speciflc plasmon modes is the least studied area of plasmonics. Nevertheless, the
temporal analysis of plasmon modes is much needed in order to fully comprehend
the time-dynamics of their excitation as well as their saturation and decay. This
temporal analysis can also be very instrumental in the area of light controllability
of plasmon resonances in semiconductor nanoparticles, where proper time synchro-
nization of excitation of speciflc plasmon modes may be required. During our future
research work on this project, the temporal analysis of plasmon modes should be
addressed.
(2) Plasmon resonances in semiconductor nanoparticles. Plasmon resonances
can be controlled through the manipulation of conduction electron density. In semi-
conductors, the manipulation of conduction electron density can be accomplished
by doping as well as by optical and depletion means. Indeed, by appropriate dop-
ing of semiconductor nanoparticles, the wide range of controllability of !p can be
achieved and, in this way, the semiconductor nanoparticles can be tuned to resonate
at desirable frequencies. The optical controllability is especially attractive because
it could be utilized for the development of nanoscale light switches and all-optical
nano-transistors.
(3) Loss reduction in plasmonic circuits. A major obstacle for the applications
of plasmonic sub-wavelength waveguides and circuits is the high losses compared
to conventional optical waveguides. The losses are mainly due to the absorption,
scattering and coupling. This problem can be attacked from two aspects: one is to
design novel metallic nanostructures to minimize the coupling losses; the other one
is to embed the nanostructures into an optically active medium to compensate for
128
the absorption and scattering losses. The successful solution of this problem could
lead to the realization of plasmonic chips.
The surface integral technique is a powerful tool for the study of plasmon res-
onances in nanoparticles and their applications. We hope that future developments
and improvement will greatly expand its capabilities and make it an indispensable
tool in this exciting and promising area.
129
BIBLIOGRAPHY
[1] National Science and Technology Council and Subcommittee on Nanoscale Sci-
ence, \National Nanotechnology Initiative: The Initiative and its Implementa-
tion Plan", Engineering and Technology Report, July 2000.
[2] S. Link and M. A. El-Sayed, \Optical Properties and Ultrafast Dynamics of
metallic Nanocrystals", Annu. Rev. Phys. Chem. vol.54, pp.331-366, 2003.
[3] U. Kreibig and M. Vollmer, Optical properties of metal clusters, Springer Series
in Material Science, vol.25, 1995.
[4] G. Mie, \A contribution to the optics of turbid media special colloidal metal
solutions", Ann. Phys. 25, pp.377-445, 1908.
[5] M. C. Daniel and D. Astruc, \Gold nanoparticles: assembly, supramolecular
chemistry, quantum-size-related properties, and applications toward biology,
catalysis and nanotechnology", Chem. Rev. vol.104, pp.293-346, 2004.
[6] B. Wiley, Y. Sun, J. Chen, H. Cang, Z.-Y. Li, X. Li, and Y. Xia, \Shape-
controlled synthesis of silver and gold nanostructures", MRS Bull., vol.30,
pp.356-361, 2005.
[7] R. C. Jin, Y. W. Cao, C. A. Mirkin, K. L. Kelly, G. C. Schatz and J. G.
Zheng, \Photoinduced conversion of silver nanospheres to nanoprisms", Sci-
ence, vol.294, pp.1901-1903, 2001.
130
[8] W. L. Barnes, A. Dereux and T. W. Ebbesen, \Surface plasmon subwavelength
optics", Nature, vol.424, pp.824-830, 2003.
[9] I. D. Mayergoyz and Z. Zhang, \Numerical Analysis of Nanoparticle-structured
Plasmon Waveguides of Light", IEEE Transactions on Magnetics, vol.43,
pp.1685-1688, 2007.
[10] S. A. Maier, P. G. Kik and H. A. Atwater, \Observation of coupled plasmon-
polariton modes in Au nanoparticle chain waveguides of difierent lengths: Esti-
mation of waveguide loss", Applied Physics Letters, vol.81, pp.1714-1716, 2002.
[11] S. A. Maier, P. G. Kik and H. A. Atwater, \Optical Pulse Propagation in Metal
Nanoparticle Chain Waveguide", Physical Review B, vol.67, 205402, 2003.
[12] S. A. Maier, P. G. Kik, H. A. Atwater, S. Meltzer, E. Harel, B. E. Koel and
Ari Requicha, \Local detection of electromagnetic energy transport below the
difiraction limit in metal nanoparticle plasmon waveguides", Nature Materials,
vol.2, pp229-232, 2003.
[13] J. C. Hulteen, D. A. Treichel, M. T. Smith, M. L. Duval, T. R. Jensen and R.
P. Van Duyne, J. Phys. Chem. B, vol.103, pp.3854-3863, 1999.
[14] C. L. Haynes and R. P. Van Duyne, \Nanosphere Lithography: A Versatile
Nano-fabrication Tool for Studies of Size-Dependent Nanoparticles Optics", J.
Phys. Chem. B, vol.105, pp.5599-5611, 2001.
131
[15] P. G. Kik, S. A. Maier, and H. A. Atwater, \Plasmon Printing A New Approach
to Near-fleld Lithography", Mater. Res. Soc. Symp. Proc., vol.705, pp.66-71,
2002.
[16] M. Moskovits, \Surface-enhanced spectroscopy", Rev. Mod. Phys., vol.57,
pp.783-826, 1985.
[17] S. Nie and S. R. Emory, \Probing Single Molecules and Single Nanoparticles by
Surface-Enhanced Raman Scattering", Science, vol.275, pp.1102-1106, 1997.
[18] K. Kneipp, Y. Wang, H. Kneipp, L. T. Perelman, I. Itzkan, R. R. Dasari,
and M. S. Feld, \Single Molecule Detection Using Surface-Enhanced Raman
Scattering (SERS)", Physical Review Letter, vol.78, 1667, 1997.
[19] D. R. Fredkin and I. D. Mayergoyz, \Resonance Behavior of Dielectric Objects
(Electrostatic Resonances)", Physical Review Letter, vol. 91, 253902, 2003.
[20] I. D. Mayergoyz, D. R. Fredkin and Z. Zhang, "Electrostatic (Plasmon) Reso-
nance in Nanoparticles", Physical Review B 72, 155412, 2005.
[21] J. Tominnaga, C. Mihalcea, D. Buchel, H. Fukuda, T. Nakano, N. Atoda, H.
Fuji and T. Kikukawa, \Local Plasmon Photonic Transistor", Applied Physics
Letters, vol.78, pp.2417-2419, 2001.
[22] R. Elghanian, J. J. Storhofi, R. C. Mucic, R. L. Letsinger and C. A. Mirkin,
\Selective Colorimetric Detection of Poly-nucleotides Based on The Distance-
Dependent Optical Properties of Gold Nanoparticles", Science, vol.277,
pp.1078-1081, 1997.
132
[23] S. Schultz, D. R. Smith, J. J. Mock, and D. A. Schultz, \Single-target mole-
cule detection with nonbleaching multicolor optical immunolabels", Proc. Natl.
Acad. Sci., vol.97, pp.996-1001, 2000.
[24] Kagan Kerman et al., \Recent Trends in Electrochemical DNA biosensor tech-
nology", Measurement Science and Technology, vol.15, pp.R1-R11, 2004.
[25] G. P. Mitchell, C. A. Mirkin and R. L. Letsinger, \Programmed assembly of
DNA functionalized quantum dots", J. Am. Chem. Soc., vol.121, pp.8122-8123,
1999.
[26] T. A. Taton, C. A. Mirkin and R. L. Letsinger, \Scanometric DNA Array
Detection with Nanoparticle Probes", Science 289, pp.1757-1760, 2000.
[27] E. Prodan, C. Radlofi, N. J. Halas and P. Nordlander, \A Hybridization
ModelforThePlasmon ResponseofComplexNano-structures", Science, vol.32,
pp.419-422, 2003.
[28] S. Lal, S. L. Westcott, R. N. Taylor, J. B. Jackson, P. Nordlander and N. J. Ha-
las, \Light Interacting Between Gold Nanoshells Plasmon Resonance and Pla-
nar Optical Wave-guides", Journal of Physical Chemistry B, vol.106, pp.5609-
5612, 2002.
[29] J. B. Pendry, A. J. Holden, W. J. Stewart and I Youngs, \Extremely Low Fre-
quency Plasmon in Metallic Meso-structures", Physical Review Letter, vol.76,
4773, 1996.
133
[30] G. Shvets and Y. A. Urzhumov, \Engineering the electromagnetic properties of
periodic nanostructures using electrostatic resonances", Physical Review Let-
ter, vol.93, 243902, 2004
[31] A. Alu, A. Salandrino, and N. Engheta, \Negative efiective permeability and
left-handed materials at optical frequencies", Optical Express, vol.14, pp.1557-
1567, 2006.
[32] J. Tominnaga, T. Nakano and N. Atoda, Proc. SPIE, vol.3467, pp.282-286,
1999.
[33] H. Ditlbacher, J. R. Krenn, B, Lamprecht, A. Leitner and F. Aussenegg, \Spec-
trally Coded Optical Data Storage by Metal Nanoparticles", Optical Letters,
vol.25, pp.563-565, 2000.
[34] L. Men, J. Tominaga, H. Fuji, Q. Chen and N. Atoda, \High-density optical
data storage using scattering-mode superresolution near-fleld structure", Proc.
SPIE, vol.4085, pp.204-207, 2001.
[35] T. J. Silva and S. Schultz \A scanning near-fleld optical microscope for the
imaging of magnetic domains in re ection", Rev. Sci. Instrum., vol.67, pp.715-
725, 1996.
[36] R. M. Stockle, Y. D. Suh, V. Deckert and R. Zenobi, \Nanoscale chemical
analysis by tip-enhanced Raman spectroscopy", Chem. Phys. Lett., vol.318,
pp.131-136, 2000.
134
[37] B. P. Rand, P. Peumans, and S.R. Forrest, \Long-range absorption enhance-
ment in organic tandem thin-fllm solar cells containing silver nanoclusters",
Journal of Applied Physics, vol.96, pp.7519-7526, 2004.
[38] N. M. Lawandy, \Localized surface plasmon singularities in amplifying media",
Applied Physics Letter, vol.85, pp.5040-5042, 2004.
[39] Craig F. Bohren and Donald R. Hufiman, Absorption and scattering of light
by small particles, John Wiley & Sons Press, 1983.
[40] I. Grigorenko, S. Haas and A. F. J. Levi, \Electromagnetic Response of Broken-
Symmetry Nanoscale Clusters", Phys. Rev. Lett., vol.97, 036806, 2006.
[41] K. L. Kelly, A. A. Lazarides and G. C. Schatz, \Computational Electromagnet-
ics of Metal Nanoparticles and Their Aggregates", Computing in Science and
Engineering, vol.3, pp.67-73, 2001.
[42] E. Hao, G. C. Schatz and J. T. Hupp, \Synthesis and optical properties of
anisotropic metal nanoparticles", Journal of Fluoresence, vol. 14, pp.331-341,
2004.
[43] E. M. Purcell and C. R. Pennypacker,\Scattering and absorption of light by
nonspherical dielectric grain", The Astrophysical Journal, vol.186, pp.705-714,
1973.
[44] B. T. Draine, \The discrete-dipole approximation and its application to inter-
stellar graphite grains", The Astrophysical Journal, vol.333, pp.848-872, 1988.
135
[45] S. K. Gray and T. Kupka, \Propagation of light in metallic nanowire arrays:
Finite-difierence time-domain studies of silver cylinders", Physical Review B
68, 045415, 2003.
[46] A Ta ove, Computational Electrodynamics: The Finite-Difierence Time-
Domain Method, Artech House, Boston, 3rd edition 2005.
[47] K. S. Yee, \Numerical solution of initial boundary value problems involv-
ing Maxwells equations in isostropic media", IEEE Trans. Antennas Propag.,
vol.14, pp.302-307, 1966.
[48] D. W. Mackowski and M. I. Mishenko,\Calculation of the T Matrix and the
Scattering Matrix for Ensembles of Spheres", J. Optical Society America A,
vol. 35, pp.2266-2278, 1996.
[49] M. I. Mishchenko, L. D. Travis, D. W. Mackowski, \T-matrix computations
of light scattering by nonspherical particles: A review", J. Quant. Spectrosc.
Radiat. Transfer, vol.55, pp.535-575, 1996.
[50] B. Peterson, S. Strom, \T-matrix formulation of electromagnetic scattering
from multilayered scatterers", Phys. Rev. D, vol.10 pp.2670-2684, 1974.
[51] E. Moreno, D. Erni, C. Hafner, and R. Vahldieck, \Multiple multipole method
with automatic multipole setting applied to the simulation of surface plasmons
in metallic nanostructures", J. Opt. Soc. Am. A, vol.19, pp.101-111, 2002.
[52] L. Novotny, D.W. Pohl, and B. Hecht,\Scanning Near-Field Optical Probe with
Ultrasmall Spot Size", Optical Letters, vol. 20, pp.970C972, 1995.
136
[53] C. Hafner, The Generalized Multipole Technique for Computational Electromag-
netics, Artech House, Boston, 1990.
[54] I. Mayergoyz, Iterative Techniques for the Analysis of Static Fields in Inhomo-
geneous, Anisotropic and Nonlinear Media, Naukova Dumka, Kyiv, 1976.
[55] O. Tozoni and I. Mayergoyz, Analysis of 3D Electromagnetic Fields, Techuika,
Kyiv, 1974.
[56] F. Ouyang and M. Isaacson, Philos. Mag. B, vol.60, pp.481, 1989.
[57] F. Ouyang and M. Isaacson, Ultramicroscopy, vol.31, pp.345, 1989.
[58] Z. Zhang and I. D. Mayergoyz, N. Gumerov, R. Duraiswami,\Numerical Analy-
sis of Plasmon Resonances Based on Fast Multipole Method", IEEE Transac-
tions on Magnetics, vol.43, pp.1465-1468, 2007.
[59] I. D. Mayergoyz and Z. Zhang,\Numerical Analysis of Plasmon Resonances
in Metallic Nanoshells with Silicon Core", IEEE Transactions on Magnetics,
vol.43, pp.1689-1692, 2007.
[60] I. D. Mayergoyz and Z. Zhang,\The Computation of Extinction Cross-sections
of Resonant Metallic Nanoparticles Subject to Optical Radiatio", IEEE Trans-
actions on Magnetics, vol.43, pp.1681-1684, 2007.
[61] I. D. Mayergoyz, Z. Zhang and G. Miano,\Analysis of Dynamics of Excitation
and Dephasing of Plasmon Resonance Modes in Nanoparticle", Physical Review
Letters, vol.98, 147401, 2007.
137
[62] I. D. Mayergoyz and Z. Zhang,\Numerical Analysis of Plasmon Resonances in
Nanoparticle", IEEE Transactions on Magnetics, vol.42, pp.759-762, 2006.
[63] I. D. Mayergoyz and Z. Zhang,\Modeling of The Electrostatic (Plasmon) Reso-
nance in Metallic and Semiconductor Nanoparticle", Journal of Computational
Electronics, vol.4, pp.139-143, 2005.
[64] L. Greengard, and V. Rokhlin, \A fast algorithm for particle simulations", J.
Comput. Phys., vol.73, pp.325-348, 1987
[65] N. A Gumerov and R. Duraiswami, Fast Multipole Methods for the Helmholtz
EquationinThreeDimensions (TheElsevierElectromagnetismSeries), Elsevier
Science, 2005.
[66] N. A. Gumerov, R. Duraiswami and E. A. Borovikov, \Data structure, op-
timal choice of parameters, and complexity results for generalized multilevel
fast multipole methods in d dimensions", UMIACS TR 2003-28, University of
Maryland, College Park, 2003.
[67] Zhihui Tang, \Fast transforms based on structured matrices with applications
to the fast multipole method", Ph.D Dissertation, University of Maryland,
2004.
[68] P. B. Johnson and R. W. Christy, \Optical Constant of The Noble Metals",
Physical Review B, vol.6, pp.4370-4379, 1972.
[69] E. D. Palik, Handbook of Optical Constant of Solids, The Academic Press, 1985.
138
[70] O. D. Kellogg, Foundation of Potential Theory, McGraw-Hill, New York, 1929.
[71] S.G.Mikhlin, Mathematical Physics, an Advanced Course, North-Holland, Am-
sterdam, 1970.
[72] J. Aizpurua and P. Hanarp and D. S. Sutherland and M. Kall and Garnett W.
Bryant and F. J. Garcia de Abajo, Phys. Rev. Lett., vol.90, 57401, 2003.
[73] K. H. Su, Q. H. Wei, X. Zhang, J. J. Mock, D. R. Smith, and S. Schultz \Inter-
particle Coupling Efiects on Plasmon Resonances of Nano-gold Particles", Nano
Letters, vol.3, pp.1087-1090, 2003.
[74] L. J. Sherry, S. H. Chang, G. C. Schatz, and R. P. Van Duyne, \Localized
surface plasmon resonance spectroscopy of single silver nanocubes", Nano Lett.,
vol.5, pp.2034-2038, 2005.
[75] W. C. Chew, J. Jin, E. Michielssen and J. Song, \Fast and E?cient Algorithms
in Computational Electromagnetics", Artech House Publishers, 2000.
[76] L. Greengard, \Fast algorithms for classical physics", Science, vol.265, pp.909-
914, 1994.
[77] V. Rokhlin, \Rapid solution of integral equations of classical potential theory",
Journal of Computational Physics, vol.60, pp.187-207, 1983.
[78] I. D. Mayergoyz, P. Andrei and M. Dimian, \Nonlinear magnetostatic calcu-
lations based on fast multipole method", IEEE Transactions on Magnetics,
vol.39, pp.1103-1106, 2003.
139
[79] L. Greengard, The Rapid Evaluation of Potential Fields in Particle System,
The MIT Press (ACM Distinguished Dissertation 1987), 1987.
[80] N. A. Gumerov and R. Duraiswami, \Multiple scattering from N spheres using
multipole reexpansion", J. Acoust. Soc. Am., vol.112, pp.2688-2701, 2002.
[81] N. A. Gumerov, R. Duraiswami and E.A. Borovikov, \Data Structures, Opti-
mal Choice of Parameters, and Complexity Results for Generalized Multilevel
Fast Multipole Methods in d Dimensions", UMIACS TR 2003-28, Also issued
as Computer Science Technical Report CS-TR-4458, University of Maryland,
College Park, 2003.
[82] N. A. Gumerov and R. Duraiswami, \Comparison of the e?ciency of translation
operators used in the fast multipole method for the 3D Laplace equation",
UMIACS TR 2005-09, University of Maryland, College Park, 2005.
[83] K. Li, M. Stockman, and D. J. Bergman, \Self-similar chain of metal
nanospheres as an e?cient nanolens", Phys. Rev. Lett., vol.91, 227402, 2003.
[84] K. Sato, K. Hosokawa, and M. Maeda, \Rapid aggregation of gold nanoparti-
cles induced by non-cross-linking DNA hybridization", J. Amer. Chem. Soc.,
vol.125, pp.8102-8104, 2003.
[85] L. Tsang, J. A. Kong, K-H. Ding, Scattering of Electromagnetic Waves |
Theories and Applications, John Wiley & Sons, 2000.
[86] G. Dassios and R.Kleinman, Low Frequency Scattering, Oxford university press,
2000.
140
[87] J. D. Jackson, Classical Electrodynamics, John Wiley and Sons Press, 3rd edi-
tion, 1999.
141