NONLINEAR SURFACE PLASMON POLARITONS: ANALYTICAL AND
NUMERICAL STUDIES
by
YAN GUO
A DISSERTATION
Presented to the Department of Physics
and the Graduate School of the University of Oregon
in partial fulfillment of the requirements
for the degree of
Doctor of Philosophy
March 2012
DISSERTATION APPROVAL PAGE
Student: Yan Guo
Title: Nonlinear Surface Plasmon Polaritons: Analytical and Numerical Studies
This dissertation has been accepted and approved in partial fulfillment of the
requirements for the Doctor of Philosophy degree in the Department of Physics
by:
Dr. Steven van Enk
Dr. Miriam Deutsch
Dr. Jens No¨ckel
Dr. Dietrich Belitz
Dr. Andrew Marcus
Chair
Advisor
Inside Member
Inside Member
Outside Member
and
Kimberly Andrews Espy Vice President for Research and
Innovation / Dean of the Graduate
School
Original approval signatures are on file with the University of Oregon Graduate
School.
Degree awarded March 2012
ii
©March 2012
Yan Guo
iii
DISSERTATION ABSTRACT
Yan Guo
Doctor of Philosophy
Department of Physics
March 2012
Title: Nonlinear Surface Plasmon Polaritons: Analytical and Numerical Studies
This dissertation contains analytical and numerical studies of nonlinear
surface-plasmon polaritons (SPPs). In our studies, we consider SPP propagation
at the interface between a noble metal with a cubic optical nonlinearity and an
optically linear dielectric.
We first consider a sum-frequency generation process during the nonlinear
interaction, where a nonlinear polarization with tripled frequency is generated from
the incident fundamental SPP. Using the non-depletion approximation, the solution
of the nonlinear wave equation shows a third harmonic generation process from the
incident SPP wave. The solution is bound in the dielectric while freely propagating
in the metal. For realistic noble metals with absorption, we use silver for its
transparency window around the plasma frequency. In this window, absorption
losses are reduced and the resultant signal has a good transmittance within the
metal. The energy conversion efficiency from the incident SPP wave to the THG
iv
signal is about 0.1% for excitation using a standard continuous wave laser with
visible light intensity I = 103W/cm2. Once generated, the propagation angle of the
signal is fully determined by the optical properties of the dielectric and the metal
layers.
We next consider a nonlinear polarization with the same frequency as
the incident light. In this process the third order nonlinearity of the metal
is described by a nonlinear optical refractive-index. With the slowly varying
amplitude approximation, the nonlinear wave equation takes the form of a
nonlinear temporal Schro¨dinger (NLS) equation. The solution to the NLS equation
for the nonlinear SPP is a temporal dark soliton (TDS). In addition to analytical
studies, computational methods are also used. With no metal loss, the numerical
solution shows stable propagation of a TDS, when the initial pulse has a tanh
envelope satisfying the threshold peak amplitude. For an arbitrary input pulse,
instabilities such as background-oscillations and multi-peak breakups occur. With
metal loss, the input optical pulse decays while maintaining a single pulse shape
when the initial amplitude satisfies the same tanh envelope condition as in the
lossless case. For an arbitrary pulse, background-oscillations or pulse-breakups
occur after a short time of propagation.
This dissertation includes previously published and unpublished co-authored
material.
v
CURRICULUM VITAE
NAME OF AUTHOR: Yan Guo
GRADUATE AND UNDERGRADUATE SCHOOLS ATTENDED:
University of Oregon, Eugene, Oregon
University of Science and Technology of China, Hefei, P.R.China
DEGREES AWARDED:
Doctor of Philosophy in Physics, 2012, University of Oregon
Master of Science in Physics, 2005, University of Oregon
Bachelor of Arts in Physics, 2003, University of Science and Technology
of China
AREAS OF SPECIAL INTEREST:
Nonlinear Optics, Nonlinear Plasmonics, Analytical and Numerical
Modeling, Numerical Computation and Simulations
PROFESSIONAL EXPERIENCE:
Graduate Research Assistant,
University of Oregon, 2005 – 2011
Graduate Teaching Fellow,
University of Oregon, 2003 – 2005
vi
GRANTS, AWARDS AND HONORS:
Travel Grant, Division of Laser Science, American Physical Society
(APS) 2010
Ph.D. Qualifying Exam First Place Award, Department of Physics,
University of Oregon 2004
PUBLICATIONS:
Yan Guo, Miriam Deutsch, “Third Harmonic Generation using Surface
Plasmon Polaritons at Nonlinear Interfaces”, Frontiers in Optics
(FiO)/Laser Science XXVI (LS) conference, 696, FWN5 (2010)
Yan Guo, Miriam Deutsch, “Nonlinear Surface Plasmon Polariton
Propagation and Third Harmonic Generation at a Planar
Metal/Dielectric Interface”, APS March Meeting Z15, (2010).
vii
ACKNOWLEDGEMENTS
My deepest gratitude is to my advisor, Dr. Miriam Deutsch, who was a
co-author of all the materials presented in this dissertation. I have been feeling
amazingly fortunate to have an advisor who gave me guidance and encouragement
during the study and research of all these years. Her patience and kindness, as well
as her enthusiasm in science and sound advice in research, have been invaluable to
me. This dissertation would not have been possible without her efforts and help.
I am deeply grateful to Dr. Peter Hugger who has provided assistance in
numerous ways. Peter has been my English guru, good roommate, and friend for
many years. Its a great pleasure to know him in and out of physics. I am also
thankful to him for his endured efforts in carefully reading and editing on countless
revisions of this dissertation.
I am indebted to many of my colleagues to support me, especially Lawrence
Mick Davis. I really cherish the countless discussions we had in research and
exchange of ideas. Its a great honor to work with him.
I wish to thank all my classmates in grad-school, Andrew, Cody, Jeremy, Matt,
Pete, Libby, Emelie and Kathy, for the study time with happiness and hardness, for
the entertainment, and for the friendship. I wish to thank my Chinese friends in the
Physics Department, Yupeng, Xianghui, Yan, Jun and Xiaolu, for making Eugene
viii
home. I also want to thank my friends in the Physics Department, Shannon, Billy,
Chuck, Mary and Carey, for helping me get through the difficult times.
This work could not have been completed without the support of the
committee: Dr. Steven van Enk, Dr. Jens No¨ckel and Dr. Dietrich Belitz in
Department of Physics, and Dr. Andy Marcus in Department of Chemistry. I
also would like to thank Dr. Susanta Sarkar for teaching me experimental skills;
Dr. Keisuke Hasegawa for the initial momentum on computer simulations and
modelings; and Bonnie Grimm and other staff members of the physics office and
the Oregon Center for Optics for backing me up inside the Willamette Hall.
Finally I would like to thank my wife, Chu Chen. She is always giving me the
strongest support and greatest understanding. Thank you for always being there
for me.
ix
To CHEN Chu and GUO Rui.
x
TABLE OF CONTENTS
Chapter Page
I. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
From Electronics and Photonics to Plasmonics . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Organization of This Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
II. INTRODUCTION TO PLANAR SURFACE PLASMON
POLARITONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Optical Properties of Noble Metals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Planar Surface Plasmon Polaritons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
III. INTRODUCTION TO NONLINEAR MEDIA AND
NONLINEAR WAVE EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
Introduction to Nonlinear Optics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
Wave Equations in a Nonlinear Medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Simplifications of the Nonlinearities in an SPP System. . . . . . . . . . . . . . . . . 38
IV. HARMONIC GENERATION IN SURFACE PLASMON
POLARITONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
Theoretical Development of Third Harmonic Generation
of Nonlinear SPPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
Conversion Efficiency and Propagation Angles . . . . . . . . . . . . . . . . . . . . . . . . . 51
Adjustment to Realistic Noble Metals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
V. PLASMONIC SOLITONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
xi
Chapter Page
Introduction to Optical Solitons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
The Temporal Nonlinear Schro¨dinger Equation of SPPs . . . . . . . . . . . . . . . 66
Nonlinear SPP Dispersion Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Analytical and Numerical Solutions of Dark Solitons. . . . . . . . . . . . . . . . . . . 77
Numerical Analysis with and without Metal Loss . . . . . . . . . . . . . . . . . . . . . . 82
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
VI. CONCLUSIONS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
A. SOLUTION OF INHOMOGENEOUS SECOND
ORDER ORDINARY DIFFERENTIAL EQUATIONS . . . . . . . . . . . . . . . . 95
B. HELMHOLTZ THEOREM AND HELMHOLTZ EQUATION . . . . . . . . . 97
C. THE ROTATION MATRIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
D. NUMERICAL ALGORITHM FOR NLS EQUATIONS . . . . . . . . . . . . . . . 113
E. C/C++ SOURCE CODE FOR FDTD METHOD . . . . . . . . . . . . . . . . . . . . 124
FDTDmain.cpp. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
FDTDarray.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
FDTDarray.cpp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
FDTDcubic.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
FDTDcubic.cpp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
FDTDtime.h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
FDTDtime.cpp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
REFERENCES CITED . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
xii
LIST OF FIGURES
Figure Page
1.1 Moore’s sketch. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Plasmonics is the bridge between electronics and photonics. . . . . . . . . . . . . . . 4
1.3 A schematic of SPPs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 All-optical switching. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 A schematic of a surface plasmon resonance biosensor. . . . . . . . . . . . . . . . . . . . 9
2.1 The modified Drude model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Sketch of an SPP-supporting interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 The dispersion relation of SPPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Prism coupling methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.5 TE and TM mode separation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.6 The penetration depths and a snapshot of the SPPs. . . . . . . . . . . . . . . . . . . . 29
3.1 Sum frequency generation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2 Difference frequency generation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3 Three coordinate systems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.1 A graph of the THG solution in vacuum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.2 A graph of the THG solution in metal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.3 A plot of the propagation angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4 Drude model vs realistic metals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.5 A plot of propagation angle with silver data. . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.1 Temporal Soliton. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2 Spatial Soliton. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3 A bright soliton and a dark soliton. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.4 The coefficients in the NLS equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.5 The stability of the soliton presented by numerical calculations. . . . . . . . . 83
5.6 Pulse break-up and oscillation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.7 The evolution of a dark SPP pulse with loss. . . . . . . . . . . . . . . . . . . . . . . . . . . 86
xiii
Figure Page
5.8 Pulse break-up and oscillation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
C.1 The rotation matrix for crystal - lab coordinates conversion. . . . . . . . . . . . 99
C.2 The relation between (111) crystal “cut” and the lab coordinates. . . . . . 102
C.3 The relation between (110) crystal “cut” and the lab coordinates. . . . . . 103
C.4 The relation between (100) crystal “cut” and the lab coordinates. . . . . . 104
D.1 Crank-Nicolson scheme illustrated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
xiv
LIST OF TABLES
Table Page
2.1 Maxwell’s equations and units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Parameters in the modified Drude model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
D.1 Numerical solution with lossless metals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
D.2 Numerical solution with lossy metals.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
xv
CHAPTER I
INTRODUCTION
From Electronics and Photonics to Plasmonics
In the computer hardware industry, the well-known Moore’s Law predicts the
growth speed of chip technology. It states that the number of transistors that can
be printed per square inch on a chip doubles approximately every 18 months, shown
in Fig. 1.1.
This famous observation was made by Intel co-founder Gordon Moore in 1965
and has described the pace of silicon technology growth ever since then. Today,
decades later, Moore’s Law remains true. However Moore’s Law is also meeting
new challenges and slowing down [1]. The major challenges opposing further
improvements are the electronic thermal effect (Joule heating) and signal speed
delay within electronic interconnections [2]. The Joule heating effect is inherent
for conductors carrying electric current, and is caused by the interactions between
the moving electrons and the background atomic ions. The heating problem shows
more dramatically when the integrated circuits enter the sub-micron dimension
in size, for both more heat is generated per unit area and the heat is not able to
dissipate effectively. The temperature built up by Joule heat in turn deteriorates the
1
Figure 1.1 Moore’s sketch. In 1965, Gordon Moore sketched out his prediction of
the pace of silicon technology. It is still true today while facing new challenges.
The picture is taken from http://download.intel.com/museum/Moores_Law/
Printed_Materials/Moores_Law_2pg.pdf. The axis labels read: “Number
of Components Per Integrated Circuit” (x-axis), and “Relative Manufacturing
Cost/Component” (y-axis)
quality and the speed of the electronic interconnections. In a world where smaller
is better, today’s state-of-the-art micro chips shrink the dimension of circuit units
to the order of 50nm [3] and the bottleneck to the development of smaller and
2
faster integrated circuit chips has become the fundamental conductivity properties
of electrons [4]. For electronic integrated circuits, the limiting factors also include
the distributed capacitance and the cross-talk between components factors which
are fundamentally caused by the finite mass of electrons, preventing the production
of electronic devices on yet smaller scales.
One of the possible solutions to the major problems in electronics is to use
photonics technology. Photons provide great data carrying capacity, and high
speed data processing photons operate at optical frequencies and travel at the
speed of light. Since direct interactions between photons are extremely weak, and
negligible in practice, photonics devices are robust to environmental photonic noise.
In addition, for a given data transmission rate, photonics devices exhibit much lower
power consumption than electronic devices and lose less energy as heat. Because
of these strengths in photonics, the invention of laser and fiber communication
technology has greatly enhanced information processing speed and accurate long
distance information relay. Presently, photonics is being investigated as a method
to connect electronic devices within sub-micron distance on a chip, as illustrated in
Fig. 1.2.
One of the fundamental limiting factors in photonics, however, is the diffraction
limit of optical waves. The typical operating wavelength in telecom optical fibers
is about 1.5µm [5], which sets the lower limit of optical devices and operation
dimension to be no smaller than hundreds of nanometers. More precisely, when the
3
Figure 1.2 Plasmonics is the bridge between electronics and photonics. It utilizes
metal/dielectric interfaces as optical wave guides and confines optical waves to
subwavelength dimensions. The operating speed of plasmonic devices is still
determined by the speed and bandwidth of light waves. Thus it takes advantage of
operating speeds of photonic devices and the small critical dimensions of electronics
devices. (Figure taken from ref. [4]).
operating wavelength of the light is λ ∼ 1500nm, the wave nature of the photons
gives the operation dimension, d, its lower bound according to Abbe diffraction
limit of light waves, d = λ/2n where n is the optical refractive index of the medium
in which the light travels. For example in glass, the operating dimension is roughly
530nm by the Rayleigh criterion since glass has a refractive index of n = 1.5. For
photons, the ability to store and transport data deteriorates when the device size
scale is under its half wavelength. Therefore photonic devices are at least one order
of magnitude larger than nano-scale electronic components.
4
In order to facilitate both optical and electronic devices integrated onto the
same chip, a bridge of interfacing needs to be built to close the gap between the
two very different technologies. This bridge is afforded by the photonics-electronics
hybrid, i.e. plasmonics (see Fig. 1.2).
The fundamental element of plasmonics is the surface plasmon-polariton
(SPP). In 1956 Pines described oscillations of the free electron gas in a conducting
metal as “plasmons”, in analogy to earlier study on plasma oscillations in gas
discharges [6]. The term “polariton” was first used by Fano in the same year to
describe the particle-like coupled oscillation formed by interactions between light
and electrons in a metal [7]. Finally the term “surface plasmon-polariton” was
coined by Cunningham and his colleagues in 1974 [8].
The interactions between light waves and the free electrons in metal lead to the
SPP modes guided by the dielectric-metal interface. As a result, the electromagnetic
wave is bound on the metal surface, with only evanescent fields into the dielectric
and metallic media. Consequently, the planar metal surface acts as a waveguide
for the surface wave and confines the electromagnetic fields to within about half
a wavelength of the surface. Even though the penetration depth in the metal is
about 10nm for SPPs, the propagation distance of the SPP is limited by a gradual
energy dissipation from the SPP into the metal, as shown in Fig. 1.3. A typical
SPP propagation length in silver or gold is on the order of microns, and can be
5
extended to millimeters with proper structural preparation of the supporting metal
[3] [9].
Figure 1.3 A schematic of SPPs. SPPs at the interface between a metal
and a dielectric material formed from electromagnetic wave and surface charge
interactions. The magnetic field is transverse only (H is in the y direction) while
the electric field has both in-plane (x and z) components. The wave can propagate
in the x direction for a distance of more than ten wavelengths before the energy is
dissipated in the metal. (Figure taken from ref. [10]).
The simplest plasmonics are constructed by an electromagnetic wave
propagating along a vacuum-metal interface as illustrated in Fig. 1.3. For more
complicated applications with plasmonics, however, more complicated interface
geometries are built for different purposes. For example, optical wave guides
using metal nanostructures have sub-wavelength critical dimension, and thus allow
further reduction of the device size to below the diffraction limit of photonics in a
variety of geometries [11] [12].
6
The scientific branch of plasmonics has been developing rapidly for decades.
Numerous applications and many variations of plasmonics design have been
explored and researched for more than fifty years. For example, a dielectric-metal-
dielectric structured stripe waveguide was designed for thermo-optic Mach-Zender
interferometric modulators and directionally-coupled switches, which utilize long-
range SPPs guided by the metal stripe [3] [13] [14] [15] [16]. In order to guide
SPPs around curvatures such as corners [12] [17] [18] without losing high efficiency
in energy transmission and switching, arrays of closely spaced metal nanoparticles
were proposed, where the arrays were used for coupled plasmon modes between the
metal particles [19]. Similarly, ultracompact modulators have been revolutionized
using plasmonics. Modulators with femtojoule switching energies and gigahertz
modulation frequencies were demonstrated using a field-effect modulation of
plasmon waveguide modes in a metal-oxide-semiconductor geometry [20] [21].
Additionally, ultrafast, broad bandwidth optical switching of SPPs has been
demonstrated. The all-optical switching of a propagating SPP signal can be
built on a metal-dielectric waveguide with direct optical modulations by optical
excitation of the metal as shown in Fig. 1.4 [22] [23] [24] .
The evanescent waves in the transverse direction of the SPPs, on the other
hand, provide sub-micron detecting sensitivities, and thus become a very important
method in gas detection, molecular optical chemical sensing and biosensoring [25]
[26] [27] [28] [29]. These sensors operate based on surface plasmon resonance (SPR)
7
Figure 1.4 A schematic showing of ultrafast all-optical switching of SPP waves
along an aluminium/silica interface. The switching is modulated by optical pump
pulses. The set up uses the grating-coupling method to excite SPPs from a laser
and to decouple signals to detectors. (Figure taken from ref. [22]).
effects. They have the great advantage of highly sensitive to the variation of the
molecule layer and fast enough for real-time detection. These sensors also have
the great advantage of being ”label-free”, as they can function without a label on
a molecule (such as radioactivity or fluorescence). One layout of such sensors is
shown in Fig. 1.5. It positions the evanescent waves of the SPP excitation on top
of a metal film on the light-incident side and the optical detector underneath. The
binding of molecules, such as O2, CO2 or DNA to the film surface changes the
optical refractive index in the immediate vicinity of the surface layer, and produces
an easily observable shift of the resonance angle of the incident light source [25] [30]
[31].
8
Figure 1.5 A schematic of a surface plasmon resonance biosensor. This biosensor
functions by monitoring the angle shift of minimum light reflection. The biosensor
operates at the critical angle that the surface plasmon resonance (SPR) condition
is satisfied and the reflection angle of the light from the surface is monitored. The
reflection angle shifts due to the change of the optical index on the flow channel
side by the amount of molecules attached at the surface. (Figure taken from ref.
[30]).
The strong confinement of the SPPs’ electromagnetic fields can greatly
enhance the nonlinearity of the host medium, and so the study of SPPs quickly
becomes also the study of nonlinear optics. Such nonlinear optical effects include
surface enhanced Raman scattering (SERS) when the localized SPP frequency is
in resonance with the incident laser frequency [32] [33]. In addition, the surface
enhanced harmonic generations, especially the second harmonic generations which
arise from the asymmetric surface profile of the SPP structure, are also observed
and demonstrated [34] [35] [36].
9
Organization of This Dissertation
This dissertation is organized as follows. In Chapter II we first review the
optical properties of the noble metals. Based on the optical propertiesf of the noble
metals, the surface plasmon polaritons (SPPs) are introduced and derived from
Maxwell’s equations. The attributes of the linear SPPs are used and expanded in
later chapters. The material in this chapter is a review of several textbooks and
previous published articles by other authors, such as [37] and [38].
Chapter III gives a general introduction to nonlinear optics, especially the
case with third order susceptibility χ(3). The third order nonlinearity has two
major branches: sum frequency generation and nonlinear optical index. Those two
branches lead to our research interest as in SPP harmonic generation (Chapter IV)
and plasmonic solitons (Chapter V). This chapter is a review of several textbooks
such as [39] [40] and [41], with original work in theoretical derivation of nonlinear
wave equations and nonlinear polarization that are specified to our case of study,
which is previously coauthored and published with Miriam Deutsch in [42].
In Chapter IV, a solution of the nonlinear wave equations is given . The
steady-state solution shows a third harmonic generation (THG) induced by the
incident SPP and can radiate through the metal. For realistic metals, experimental
data adapted from the literature are used. The realistic metal is chosen to be silver
because silver has a transparency window around its plasma frequency and gives
good transmittance for the induced THG wave. The energy conversion efficiency
10
from the incident SPP to the THG is calculated and shown to be consistent with the
assumptions and approximations made for the original nonlinear wave equations.
The propagation angle for the THG is also derived, using both an analytical model
of the metals and experimental data from the literature. The material in this
chapter was previously coauthored and published with Miriam Deutsch in [42] and
[43].
In Chapter V, we follow established approaches for our nonlinear SPP system
and achieve a temporal nonlinear schro¨dinger (NLS) equation. The solution to this
NLS equation gives a temporal dark soliton. Attributes such as stability is studied
under scenarios of both lossless and lossy metal, with standard numerical methods
and programming. The material in this chapter was previously coauthored with
Miriam Deutsch and has not yet unpublished
11
CHAPTER II
INTRODUCTION TO PLANAR SURFACE PLASMON POLARITONS
Optical Properties of Noble Metals
The physics of plasmonics derives from electromagnetic modes called surface
plasmon polaritons (SPPs), which arise from the interaction between light and
mobile surface charges in the metal. The SPPs usually propagate along a
conductor-dielectric interface and the conductors are commonly metals with high
conductivities and low absorption such as Cu, Al and some noble metals (Ag and
Au) [10] [44] [45].
Before we dive into the world of surface plasmon polaritons, the optical model
of the main supporting medium, the noble metals, needs to be discussed in detail.
In the field of classic electromagnetism, Maxwell’s equations are the cornerstone
of all phenomena and analysis [38]. We start our analysis with the macroscopic
Maxwell’s equations in the following form:
∇ ·D = ρfree (II.1a)
∇ ·B = 0 (II.1b)
∇× E = −∂B
∂t
(II.1c)
∇×H = ∂E
∂t
+ Jtotal (II.1d)
12
We use SI units, see Tab.2.1.
Table 2.1 Maxwell’s equations and units
Symbol Name Unit
E Electric Field V/m
B Magnetic Field T
D Electric Displacement Field C/m2
H Auxiliary Magnetic Field A/m
ε0 Vacuum Permittivity F/m
µ0 Vacuum Permeability N/A
2
For the current density Jtotal, the most general expression [39] in terms of
multipole expansion is
Jtotal = Jc +
∂P
∂t
+∇×M− ∂
∂t
(∇ ·Q) + · · · (II.2)
where P is the electric dipole, M is the magnetic dipole, Q is the electric
quadrupole, and so on. In the absence of sources, ρfree = 0 the first term in (II.2),
the conducting current density, disappears, Jc = 0. In non-magnetic media, which
react only slightly to magnetic fields, the magnetization M in the third term of
(II.2) is considered negligible in optical interactions. Consequently, the relative
permeability of non-magnetic media, for our case the noble metals, is very close to
one in B = µ0µH and we can safely use µ = 1 in the context of optical properties
of the noble metals. This leads to a simple linear relation between the magnetic
induction and magnetic field: B = µ0H. For terms related to the electric charge
distribution in (II.2), the dipole term
∂P
∂t
dominates and the rest of the higher
13
order terms - quadrupole Q and higher - are dropped. Thus the optical properties
in a medium are determined to the first order of the electric field E by the dipole
polarization P and are also represented in the dielectric displacement D as
D = ε0E+P (II.3)
The relationship between the dipole polarization P and the electric field E is
macroscopically determined by the dielectric susceptibility, χ, via
P = ε0χE (II.4)
where χ is a tensor but is often simplified to a scalar quantity for isotropic media,
such as noble metals.
In a certain frequency limit, usually around the ultraviolet frequency range,
the optical response of the metal permittivity can be described by the plasma
model, where a free electron gas with number density n¯ and charge e moves against
positive background made of ion cores. In the plasma model, the electron mass me
is corrected to an effective optical mass m due to the band structure of the medium
[46]. The fact that most metals possess very high electric conductivities is a direct
result of the plasma model and we can use the plasma model for most of our studies.
In the plasma model, the conduction electrons in the metal can respond and follow
the frequency of the external electric field up to a limit, and this cut-off frequency
ωp =
√
n¯e2
mε0
is known as the plasma frequency. For alkali metals where the single
conduction electron per primitive cell is simply treated as completely free, the
14
plasma frequency ωp occurs above ultraviolet optical frequencies. For noble metals,
on the other hand, due to the presence of d−bands, the plasma frequency is in the
visible frequency range, where absorption from interband transitions dominates
[46]. We take this effect into account in later chapters. For the purposes of this
research, silver and gold are the best candidates to support SPPs. Silver and gold
are members of the noble metals family. They are chosen for both their chemical
stability and electric conductivity. Henceforth the “noble metals” referred to in
this document are silver or gold unless otherwise stated.
The electric permittivity, or dielectric function, is the primary determining
factor of optical properties of the noble metals. Here we follow the classic electron
gas model for the dielectric function [47], known as the Drude-Lorentz model for
dielectrics. The electrons in silver and gold take the effective mass, m = 1.1me
where me is the mass of an electron. They obey Newton’s second law of motion in
the presence of an external electric field E(x, t):
−eE = m[x¨ + γpx˙+ ω20x] (II.5)
In the Drude model for metals, the restoring force is considered zero for
mobile electrons such that ω0 = 0. The motion of the free electrons in the noble
metals is damped via collisions with crystal impurities and is characterized by
a phenomenological damping constant γp = 1/τrelax ∼ 100 THz where τrelax is
the relaxation time of the free electron gas on the order of femtosecond at room
temperature [46].
15
A harmonic solution x(t) = x0e
−iωt can be written down for (II.5) with an
electro-magnetic field of the form E(x, t) = E(x)e−iωt, where x0 is a complex
amplitude representing the phase shift between x and E. With the solution
x(t) =
e
m(ω2 + iγpω)
E(t) (II.6)
it is straightforward to show the macroscopic polarization has the form
P = −Nex = − Ne
2
m(ω2 + iγpω)
E (II.7)
where n¯ is the electron density in metal. Using the expression for electric
displacement (II.3) yields
D = ε0
(
1− ω
2
p
ω2 + iγpω
)
E (II.8)
where ω2p ≡
Ne2
ε0m
. So the dielectric function for a noble metal is
ε(ω) = 1− ω
2
p
ω2 + iγpω
(II.9)
The positive background of the ion cores in the metal induces a residual polarization
for ω → ∞, therefore a dielectric constant ε∞ is introduced to the Drude-Lorentz
model and the dielectric function for metal is written as
εm(ω) = ε∞ −
ω2p
ω2 + iγpω
. (II.10)
Empirically, gold and silver exhibit much stronger optical absorption than the
Drude model predicts due to the effect of the actual band structures in the metals
[48] [49] and the skin depth change in the measured experiment data [50] [51].
16
Consequently an effective absorption term i
σ
ε0ω
is added to better fit the model
where σ is the effective AC conductivity at optical frequencies [52] [53]. In this
dissertation, when absorption is important to the results, the modified Drude model
is used
ε(ω) = ε∞ −
ω2p
ω2 + iγpω
+ i
σ
ε0ω
. (II.11)
For silver, we mostly follow [54] and the parameters are listed in Table[2.2].
Table 2.2 Parameters in Modified Drude Model [37].
Symbol Value Unit
ε∞ 5.1 −
ωp 8.6 eV
γp 45 meV
σ 23.35 (Ω · µm)−1
The Drude model represents the optical properties of silver very well below the
plasma frequency ωp, whose wavelength equivalent is about 320nm for wavelength
equivalence, but fails for frequencies above ωp. A comparison of the real and
imaginary parts of ε taken from both data and the theoretically predicted values
given by the modified Drude model is given in Fig. 2.1. The experimental data and
theoretical results from the Drude model agree well over all wavelengths except
those approaching the plasma frequency. The difference between the experimental
data and the Drude model will be discussed in Ch.IV.
17
0 200 400 600 800 1000 1200
-50
-40
-30
-20
-10
0
10
 
 
R
e
[
(
)
]
Wavelength (nm)
0 200 400 600 800 1000 1200
0
1
2
3
4
5
 
 
I
m
[
(
)
]
Wavelength (nm)
Figure 2.1 The modified Drude model (lines, from (II.11)) vs experimental data
(dots, taken from [37]) for silver. For wavelengths greater than 320nm, the
Drude model correctly represents the optical properties of silver. For wavelengths
shorter than 320nm, however, the increasing absorption due to inter-band transition
dominates and the Drude model is invalid.
Planar Surface Plasmon Polaritons
Introduction
The supportive structure of surface plasmon polaritons requires at least one
interface formed by a metal and a dielectric, as shown in Fig. 2.2. The y direction
is extended to infinity in our model so that the planar SPPs discussed here are
essentially two-dimensional, with only x and z dependence.
The SPPs, which are EM waves that are guided by a metal surface, have
a larger propagation wavevector than the wavevector in vacuum for the same
frequency. Therefore additional momentum is required to excite an SPP wave
18
xz y
Dielectric
Metal
Figure 2.2 Sketch of an SPP-supporting interface. x is the direction of
propagation, z is the normal to the interface pointing from the metal to the
dielectric, and y is the transverse direction extended to infinity.
along a metal surface as seen in Fig. 2.3. There are two classical experimental
configurations used to excite the SPP from light in free space, the Kretschmann
configuration and the Otto configuration (see Fig. 2.4). Both configurations utilize
the fact that the wavevector of light is larger in a high index prism than in free
space, so when the total internal reflection condition is satisfied, the evanescent
wave tunnels out and excites SPPs on the metal surface.
There are also ways to excite SPPs without a high index prism, such as by
scattering light over a rough surface, or by using manufactured gratings to match
the missing momentum between light in free space and the SPP on a planar metal
surface, as shown in Fig. 1.4. Recently it was also demonstrated that one may excite
the SPP directly from free space by optical nonlinearities in the SPP supporting
dielectric or metal [55] [56].
The simplest geometry layout that supports surface plasmons is shown in
Fig. 2.2. The surface plasmon polaritons supported by this layout are therefore
19
kω
vacuum
ω = ck silica
ω = ck/n
SPP
∆k
Figure 2.3 The dispersion relation of SPPs. There is a momentum mismatch (red
arrow) between the SPP wavevector (blue solid curve) and the vacuum light line
(black solid line). Therefore, some special techniques are needed to compensate for
the missing momentum, such as using a high index prism (green dashed line) to
alter the light line to intersect (red dot) with the SPP dispersion curve.
called planar SPPs, differentiating from SPPs supported by other geometrical
structures such as nanoshells [57]. In photonics and plasmonics, the planar surface
layout is frequently used in chemical and biological applications such as biosensing
and nanoscale waveguides [27] [29] [30] [25].
Derivation of SPPs Using Maxwell’s Equations
We start the discussion with the Maxwell equations:
Combining the two curl equations (II.1c) and (II.1d) gives
∇×∇×E = −µ0∂
2D
∂t2
(II.12)
20
Kretschmann Otto
Figure 2.4 Prism coupling methods. Prism coupling to SPPs using attenuated
total internal reflection. The Kretschmann configuration (left) is arranged as:
prism-metal-dielectric and the Otto configuration (right) is arranged as: prism-
dielectric-metal.
in which we use the fact that the spatial and temporal derivatives are interchangeable.
Next, for linear media, optical properties can be described by the relative
permittivity ε giving D = εε0E. Considering the fact that ε0µ0 = 1/c
2, then
(II.12) is
∇×∇× E = − ε
c2
∂2E
∂t2
(II.13)
Since we have 0 = ∇ ·D = ε0∇ · (εE) = ε0E · ∇ε + ε0ε∇ · E, where the term ∇ε
is negligible when the spatial variance of ε is small over one optical wavelength, we
have in general
∇2E− ε
c2
∂2E
∂t2
= 0 (II.14)
Since the time dependence of the eigensolution is e−iωt, we write E(x, y, z, t) =
E(x, y, z)e−iωt to obtain the Helmholtz equation
(∇2 + k20ε)E(x, y, z) = 0 (II.15)
21
where k0 ≡ ω/c is the wavevector of the light with frequency ω in vacuum. Next,
as depicted in Fig. 2.2, x is the direction of propagation, z points in the direction
normal to the metal surface and y is the metal surface transverse extension such
that E has no dependence on y, i.e. ∂yE = 0. The electric field can be written as
E(x, y, z) = E(z)eiβx where β is the propagation constant. Along the transverse
direction, the form of the eigensolution satisfies ∂zE(z) = −ikzE(z). Since ∇2 =
∂2x + ∂
2
y + ∂
2
z , so it holds that
k2z + β
2 = k20ε (II.16)
Here we summarize the eigenvalues of spatial derivatives:
∇2E = −k20εE
∂xE = iβE
∂yE = 0
∂zE = ikzE
Now the dielectric function ε = ε(z, ω) is written explicitly as a piecewise constant
function with a discontinuity at the dielectric-metal interface z = 0,
ε(z, ω) =

εd(ω) when z > 0
εm(ω) when z < 0
(II.17)
where subscripts d and m denote dielectric and metal respectively. The wavevector
along the direction of propagation must be continuous across the interface, so β
22
remains the same for both media. The z component of the wavevector, however,
differs in the two media
kz(z) =

kdz when z > 0
kmz when z < 0
(II.18)
The electric field can now be written in the form
E(x, y, z, t) = E0e
iβxeikzze−iωt (II.19)
The physical attributes of the SPPs are clearly represented by the wavevectors (β
and kz), and the wavevectors can be determined in terms of the EM wave frequency
ω via the medium’s optical response - the dielectric function. Since the dielectric
function ε(ω) is a function of wave frequency too, essentially we are looking for
the dispersion relation of the eigensolution in the form of k = k(ω). To get the
SPP dispersion via a relatively straightforward path, we step back to the general
source-free Maxwell equations and then consider the boundary conditions with a
specific medium structure.
The Maxwell equations in vector form are
∇ ·D = 0 (II.20a)
∇ ·B = 0 (II.20b)
∇× E = −∂B
∂t
(II.20c)
∇×H = ∂D
∂t
(II.20d)
23
The first two divergence equations provide the boundary conditions which will be
used shortly and the last two curl equations provide the relations between different
field components. To implement the SPP geometry layout, the last two equations
can be rewritten separately by their Cartesian components in frequency ω domain
iωBx = ∂yEz − ∂zEy (II.21a)
iωBy = ∂zEx − ∂xEz (II.21b)
iωBz = ∂xEy − ∂yEx (II.21c)
−iωDx = ∂yHz − ∂zHy (II.21d)
−iωDy = ∂zHx − ∂xHz (II.21e)
−iωDz = ∂xHy − ∂yHx (II.21f)
In order to separate the 12 field components into two separate mode groups: we
distinguish transverse electric (TE) modes and transverse magnetic (TM) modes.
The SPP supporting media have linear susceptibilities with scalar response
so that D ‖ E and H ‖ B. Further, the linearity dictates that the Cartesian
components in the electric displacement D field are solely determined by the
corresponding electric E field components, as indicated by the out-going arrows
(black) in Fig. 2.5. Similarly the Cartesian components in the H field are solely
determined by the corresponding B components, as indicated by the in-going arrow
(black) in Fig. 2.5. Formally speaking, this linear dependency between the field
pairs D − E and H − B need not be assumed. In such instances, the dielectric
24
Dx
DzDy
Ex
EzEy
TE
Hx
HyHz
Bx
By
TM
Bz
∂yE
=
0
∇×E = −∂B
∂t
(outer hexagon)
∇×H = ∂D
∂t
(inner hexagon)
Figure 2.5 TE and TM mode separation. The outer hexagon shows that every B
component is determined by two E components (blue arrows) and the inner hexagon
shows that every D component is determined by two H components (pink arrows).
The out-going and in-going arrows (black) indicate linear relations in D − E and
H−B field pairs. The independence of all the fields in the y direction separate the
field components into two mode sets: the Transverse Electric (TE) mode and the
Transverse Magnetic (TM) mode with x the direction of propagation.
permittivity takes a diagonal tensor form. For this study, however, the chosen
dielectric and the metal are isotropic media, so a scalar ε suffices.
For our SPP geometry choice, a geometrical constraint of the system is that the
derivatives of the fields respect to y are zero. Since every component of the magnetic
field, B, is determined by two components in the electric field E, and similarly every
component of the electric displacement field D is determined by two components
25
of the auxiliary magnetic field H, via partial derivatives, as indicated by the outer
arrows (blue) and inner arrows (pink), the condition ∂yE = 0 essentially detaches
the x ↔ z dependency in E → B and H → D as depicted by the dashed line in
Fig. 2.5.
The eigenmode involving fields By, Ex and Ez is called Transverse Magnetic
(TM) because the magnetic field By (or Hy) is in the transverse direction relative
to the propagation direction x, and the eigenmode involving fields Ey, Hx and Hz
is called Transverse Electric (TE) for Ey transverse to x:
TM mode:(Hy, Ex, Ez)
TE mode:(Ey, Hx, Hz)
The general boundary conditions from Maxwell’s equations are
n · (D2 −D1) = ρe (II.22a)
n · (B2 −B1) = 0 (II.22b)
n× (E2 − E1) = 0 (II.22c)
n× (H2 −H1) = Ke (II.22d)
where n = zˆ is the normal direction on the surface pointing from the metal into
the dielectric.
Again, for the SPP configuration we consider, the surface charge density ρe
and the surface current density K are zero. The boundary conditions in terms of
26
Cartesian components at boundary z = 0 are:
Ddz −Dmz = 0 (II.23a)
Bdz − Bmz = 0 (II.23b)

Edx −Emx = 0
Edy − Emy = 0
(II.23c)

Hdx −Hmx = 0
Hdy −Hmy = 0
(II.23d)
here all fields marked with a subscript “d” are evaluated at z = 0+ (the vacuum
side) and all fields marked with subscript “m” are evaluated at z = 0− (the metal
side).
Explicitly using (II.20) and (II.23) we derive that:
Bx : kdzEy = kmzEy (II.24a)
By : kdzEx − βEdz = kmzEx − βEmz (II.24b)
Bz : βEy = βEy (II.24c)
Dx :
kdzHy
εd
=
kmzHy
εm
(II.24d)
Dy :
kdzHx − βHz
εd
=
kmzHx − βHz
εm
(II.24e)
Dz : βHy = βHy (II.24f)
To satisfy (II.24) requires Ey = 0 (for example (II.24a)), so the eigenmode can only
27
be TM with Hz, Ex, Ey non-zero. The relation turns out:
kdz
εd
=
kmz
εm
(II.25)
With the relation (II.16) β2 + k2z = k
2
0ε, the dispersion relation is derived:
β2 = k20
εdεm
εd + εm
(II.26a)
k2dz = k
2
0
ε2d
εd + εm
(II.26b)
k2mz = k
2
0
ε2m
εd + εm
(II.26c)
We also note the following relations between the field components from (II.24b):
Edz = −kmz
β
Ex (II.27a)
Emz = −kdz
β
Ex (II.27b)
Using the relations above, we see that the SPP profile is completely determined
by its frequency (or wavelength) and the optical response of the supporting
dielectric-metal.
In the ideal case, the metal is non-absorptive and so the dielectric function is
a real negative number. Thus the propagation factor β is real and positive while
the transverse decay factors kdz and kmz are purely imaginary. In reality, due to
the finite absorption in a metal, the dielectric function of the metal has a non-zero
imaginary component, therefore all quantities are complex numbers. For a laser
source with wavelength λ = 800nm, the corresponding SPP propagation length
Lspp and penetration lengths in the dielectric δd, and in the metal δm can be easily
28
derived from the expressions of SPP wavevectors derived above. This is shown
schematically in Fig. 2.6a for the following values:
Lspp = 1/Im[β] = 80µm
δd = 1/Im[kdz] = 635nm
δm = 1/Im[kmz] = 25nm
A real-time snapshot of the SPP field profile is shown in Fig. 2.6b, calculated
using a Finite-difference time-domain (FDTD) simulation.
Figure 2.6 The penetration depths and a snapshot of the SPPs. (a) Penetration
depth of an SPP in a dielectric and metal (Figure taken from ref. [10]) (b) The
SPP instantaneous profile snapshot, represented by the magnetic H = yˆHy field
only. This picture is achieved by an FDTD simulation written by the author, and
the source code is attached in Appendix E.
29
CHAPTER III
INTRODUCTION TO NONLINEAR MEDIA AND NONLINEAR WAVE
EQUATIONS
Introduction to Nonlinear Optics
This work was published in APS March Meeting in 2010. Dr. Miriam Deutsch
set up the proposed frame of the work I was the primary contributor to all the
derivations and analysis of the results.
Before the invention of the laser, light sources used in laboratories were
of relatively low intensity, therefore most optical phenomena observed in labs
were linear, and most experimental results could be explained simply by the
superposition principles. In reality, however, The optical properties of a material
system can be modified by the presence of high intensity light. This light-induced
modification was first demonstrated in the second harmonic generation experiment
by Franken et al. in 1961 [58]. This experiment, which occurred only one year
after the demonstration of the first working laser in 1960, marked the birth of a
new field of research in physics: nonlinear optics.
A lot of materials can show nonlinearity at different thresholds of light
intensity. Even in vacuum photons can interact via vacuum polarizations and show
30
a certain level of nonlinearity in the context of wave mixing [59]. In practice, certain
dielectric media are preferred such as glass and semiconductors [60] to observe
nonlinear optical effects. This is because nonlinear effects are greatly enhanced by
the medium polarization. Numerous new nonlinear optical phenomena have been
discovered in vast types and states of materials [61] [62] [63] [64]. The microscopic
origin of the nonlinearities in materials is classically illustrated in the nonlinear
oscillator model [39] [40].
Metals, contrary to our intuition, possess very high values of intrinsic
nonlinearity even compared to fiber glasses (SiO2) [41]. It is rare to notice the
metallic nonlinearities, however, due to the limited penetration depth of light into
the bulk metals which is typically less than one tenth of a wavelength of the light
wave. Therefore the nonlinear interactions between the light and metals are very
weak because the volume of material available to interact with the light is small.
The SPP waves, on the other hand, can concentrate high intensity of light onto the
metal surface. The SPPs also greatly increase the interaction volume between the
light and the metal by the guided nature of the SPP waves. Therefore, the SPPs
are one of the best candidates to observe nonlinearity of the metals. There have
been already different orders of metal nonlinearities demonstrated by harmonic
generations [65] [66] and by four wave mixing [56] [67].
31
Wave Equations in a Nonlinear Medium
The nonlinearity of medium can be directly expressed by means of the
macroscopic polarization P, which appears in the definition of electric displacement
D, as shown in (II.3). Under the dipole approximation, which is the lowest order
expression of P which neglects the magnetic terms, the polarization P can be
expressed as an expansion of the electric field E in polynomial orders. The
relations between polarization P(ω), electric displacement D and electric field
E(ω) are
P(ω) = ε0(χ
(1)E+ χ(2)EE+ χ(3)EEE+ · · · ) (III.1a)
D(ω) = ε0E+P = ε0(E+ χ
(1)E+ χ(2)EE+ χ(3)EEE+ · · · ) (III.1b)
where the χ(n)(ω) terms are the susceptibilities of different orders. The first two
terms, E and χ(1)E are recognized as the linear parts and can be merged and
expressed by the dielectric function:
ε(ω) = 1 + χ(1)(ω) (III.2)
and the rest is named as nonlinear polarization:
PNL(ω) = ε0(χ
(2)EE+ χ(3)EEE+ · · · ) (III.3)
Gather the equations above, together with (II.12), and the nonlinear wave equation
in a nonlinear medium is achieved with the form:
∇×∇× E+ ε
c2
∂2E
∂t2
= −µ0∂
2PNL
∂t2
(III.4)
32
In general, PNL = PNL(E) is a function of electric fieldE expanded with polynomial
terms. In order to simplify and solve the nonlinear wave equation, the detailed
nonlinear properties of the material in which the light propagates needs to be
discussed and expressed in a specific form of a response function.
For noble metals such as silver and gold which possess inversion symmetry,
the even orders of χ(n) vanish and the dominating nonlinear term is then χ(3) [39].
The nonlinear polarization is now simplified to
PNL(ω) = ε0χ
(3)(ω : ωα, ωβ, ωγ)E(ωα)E(ωβ)E(ωγ) (III.5)
where subscripts α, β and γ denote the input electric field distinguishable by their
angular frequencies ω or input momentum vectors k, and the resultant polarization
frequency is a linear combination of ωα, ωβ and ωγ, for the energy is conserved.
This is generally known as a parametric process.
In our study, we use degenerate input with the same frequency, i.e. ωa = ωb =
ωc = ω. Among several possibilities, there are two possible resultant frequencies:
3ω(= ω + ω + ω) Fig. 3.1 and ω(= ω + ω − ω) Fig. 3.2. The first case with 3ω is
introduced as tripled frequency in section .1 and discussed in details in Chapter IV
as third harmonic generation. The second case with ω is introduced as nonlinear
optical index in section .2 and the details of the resultant soliton solution are
discussed in Chapter V.
Sum Frequency - Harmonic Generation
One of the resultant nonlinear polarizations has a frequency that comprises
33
ωα
ωβ
ωγ
ω = ωα + ωβ + ωγ
χ(3)
(a)
ωα
ωβ
ωγ
ω
(b)
Figure 3.1 Sum frequency generation. (a) Sum frequency generation via nonlinear
material. Three incident beams with frequencies ωα, ωβ and ωγ interact within the
nonlinear medium with third order susceptibility χ(3). A resultant sum frequency
signal ω = ωα + ωβ + ωγ is generated.(b) Schematic illustration of energy levels for
the sum frequency process. Three photons at lower frequency are destroyed and a
photon is generated at a the frequency of sum of the destroyed photons.
the summation of the frequencies of the input light waves, known as the sum
frequency process, where three photons are converted into a new photon with much
higher energy (or frequency) Fig. 3.1b. For a sum frequency process with the same
degenerate frequency 3ω = ω + ω + ω, the nonlinear polarization term is now:
Pi(3ω) = ε
∑
jkl
χ
(3)
ijkl(3ω)Ej(ω)Ek(ω)El(ω) (III.6)
34
ωα
ωβ
ωγ
ω = ωα + ωβ − ωγ
χ(3)
(a)
ωα
ωβ ωγ
ω
(b)
Figure 3.2 Difference frequency generation. (a) Difference frequency generation
via nonlinear material. Three incident beams with frequencies ωα, ωβ and ωγ
interact within the nonlinear medium with third order susceptibility χ(3), a resultant
difference frequency signal ω = ωα+ωβ−ωγ is generated.(b) Schematic illustration
of energy levels for the difference frequency process.
and this special case of the general sum frequency processes is also called third
harmonic generation.
Now the electric field E in (III.4) has two frequency components and may be
written in the form:
E = E(ω)e−iωt + E(3ω)e−3iωt (III.7)
Therefore the nonlinear equation (III.4) is Fourier transformed into the frequency
domain by t→ ω:
∇×∇× [E(ω)e−iωt + E(3ω)e−3iωt]− ε
c2
[ω2E(ω)e−iωt + (3ω)2E(3ω)e−3iωt]
= µ0[ω
2PNL(ω)e−iωt + (3ω)2PNL(3ω)e−3iωt]
(III.8)
It is obvious that the electric and polarization fields dependent on ω and
35
3ω are coupled via (III.6) and (III.8), but the two equations can not be solved
simultaneously with exact analytical methods. In order to further simplify the
equations, the fields with the same frequency dependence are grouped together and
separated from other frequency components. Therefore the non-depletion (ND)
approximation is used to decouple the equations with different frequencies. The
ND approximation states that the nonlinear polarization drains very little energy
from the input source and that essentially the change of the input field is negligible.
The ND approximation is valid only for high input with low efficiency output. The
ND approximation insures the nonlinear polarization is solely determined by (III.6)
and will not affect the input electric field E(ω). This essentially eliminates the
dependence of between the 3ω and ω components of the fields in (III.6). Now (III.6)
can be written as two separate equations by matching the oscillation frequencies.
The one with 3ω frequency is
∇×∇×E(3ω)− ε
c2
(3ω)2E(3ω) = µ0(3ω)
2PNL(3ω) (III.9)
After the nonlinear polarization PNL(3ω) is determined from the input fields
E(ω), the generated harmonic fields E(3ω) are solved from (III.9). The incident
linear SPP wave preserves its longitudinal and transverse profiles within the ND
approximation, thus we are looking for a steady-state solution for the output field
in the form E(r, t) = E(r)eωst with a steady-state solution frequency ωs = 3ω0.
We ignore the transient solutions. The validity and standard procedure of the
steady-state solution approach are discussed in Appendix A.
36
Difference Frequency - Nonlinear Optical Index
The optical refractive index of many materials depends on the intensity of the
incident light and can described by the relation
n = n0 + n2|E|2 (III.10)
where n0 is the linear refractive index and n2 is the nonlinear index of refraction
[41]. The phenomena of a light intensity dependent optical index is sometimes
referred to as the optical Kerr effect. The nonlinear polarization in the nonlinear
Kerr effectis determined by
PNL(ω) = 3χ(3)(ω = ω + ω − ω)|E(ω)|2E(ω) (III.11)
where the intensity I ∼ |E|2.
The expression of the nonlinear polarization is a special case of (III.5) in
difference frequency generation Fig. 3.2 when the frequencies of the incident light
degenerate ωα = ωβ = ωγ = ω. Since the resultant nonlinear polarization has
the same frequency of the incident electric fields, the ND approximation for the
frequency decoupling method from last section is not applicable, no matter how
small the nonlinear polarization. The induced and original electric fields are
indistinguishable and the original beam is treated as it undergoes a nonlinear
refractive index. A new approach is needed to solve the nonlinear wave equation
and we turn to some other approximations and techniques, known as the Nonlinear
37
Scho¨dinger (NLS) Equation, which will be introduced and discussed in detail in
Chapter V.
Simplifications of the Nonlinearities in an SPP System
In the previous sections, the temporal property of the third order susceptibility
χ(3) was separated into two different frequency cases - sum and difference frequency
processes. In this section we discuss the spatial attributes of the third order
susceptibility χ(3). This discussion will be based on the sum frequency process, but
the same theory can be applied to the difference frequency process straightforwardly.
The third order susceptibility χ(3) has a general form of a rank 4 tensor with 81
components. For example, the sum frequency polarization component in Cartesian
coordinates is written:
Pi(ω : ωa, ωb, ωc) = ε0
∑
jkl
∑
(abc)
χ
(3)
ijkl(ω : ωa, ωb, ωc)Ej(ωa)Ek(ωb)El(ωc)
= ε0g
∑
jkl
χ
(3)
ijkl(ω)E
a
jE
b
kE
c
l (III.12)
The subscripts i, j, k, l represent the permutations of Cartesian axis x, y and z,
with the distinguishable input fields labeled with a, b and c. The g factor is
the degeneracy factor for all different frequency combinations and contains all the
information about frequency mixing discussed in previous sections. For example in
the degenerate case ωα = ωβ = ωγ = ω, g(3ω) = 1 since there is only one way to
38
produce the third harmonic generation
3ω = ω + ω + ω (III.13)
while g(ω) = 3 since there are three ways to produce the same frequency of ω:
ω = ω + ω − ω
ω = ω − ω + ω
ω = −ω + ω + ω
(III.14)
There ar many independent tensor elements in χ(3). Fortunately, many of
these tensor elements can be eliminated by symmetry and proper choice of the
coordinate system. Since the crystal structures of gold and silver are both face
centered cubic (FCC), the independent elements in the tensor is reduced to 21
nonzero elements with only 4 of them independent [68] [69]: χ1111, χ1122, χ1212, χ1221,
where the subscripts 1 and 2 stand for any two orthogonal Cartesian axis parallel to
the optical axis of the crystal lattice, for example χxxyy = χyyzz = χzzyy = χzzxx =
· · · = χ1122.
There are three different coordinate systems in which the electric field is
represented: the coordinates attached to the optical axis of the cubic crystal lattice
(cubic, labeled by Greek symbols α, β and γ), the orientation of the crystal lattice
(crystal, labeled by letters with a prime ′), and the lab coordinates (lab, labeled
by letters) where actual experimental measurement occurs. A simple illustration of
the coordinate systems in 2-dimension is presented in Fig. 3.3.
39
α (cubic)
γ (cubic)
x′ (crystal)
z′ (crystal)
E
Eα
Eγ
x (lab)
z (lab)
Figure 3.3 Three coordinate systems are illustrated here. The components of the
electric field E are represented relative the optical axis (cubic). For each coordinate
system, only two axes are shown for clarity.
The nonlinear polarization PNL can be written in component form, in both
cubic and lab Cartesian coordinate systems:
PNL = αˆPα + βˆPβ + γˆPγ = xˆPx + yˆPy + zˆPz (III.15)
The nonlinear polarization PNL and the electric field E have straightforward
relations in the cubic coordinate system:
PNLα = χ1111E
a
αE
b
αE
c
α
+ χ1122E
a
α(E
b · Ec)− χ1122EaαEbαEcα
+ χ1212E
b
α(E
c ·Ea)− χ1212EaαEbαEcα
+ χ1221E
c
α(E
a ·Eb)− χ1221EaαEbαEcα
(III.16a)
40
PNLβ = χ1111E
a
βE
b
βE
c
β
+ χ1122E
a
β(E
b ·Ec)− χ1122EaβEbβEcβ
+ χ1212E
b
β(E
c · Ea)− χ1212EaβEbβEcβ
+ χ1221E
c
β(E
a ·Eb)− χ1221EaβEbβEcβ
(III.16b)
PNLγ = χ1111E
a
γE
b
γE
c
γ
+ χ1122E
a
γ (E
b · Ec)− χ1122EaγEbγEcγ
+ χ1212E
b
γ(E
c · Ea)− χ1212EaγEbγEcγ
+ χ1221E
c
γ(E
a · Eb)− χ1221EaγEbγEcγ
(III.16c)
Together, in cubic coordinates (in terms of α, β, γ)
PNL = χ1122E
a(Eb · Ec) + χ1212Eb(Ec ·Ea) + χ1221Ec(Ea ·Eb)
+ χd(αˆE
a
αE
b
αE
c
α + βˆE
a
βE
b
βE
c
β + γˆE
a
γE
b
γE
c
γ)
(III.17)
where the first three terms are invariant under crystal rotations, while the last term
with χd ≡ χ1111 − χ1122 − χ1212 − χ1221 shows anisotropy under crystal rotations.
Expressing the polarization in a compact form under lab coordinate yields:
PNL = χ1122(E
b ·Ec)Ea + χ1212(Ec · Ea)Eb + χ1221(Ea · Eb)Ec + χdM :: EaEbEc
(III.18)
Therefore the tensor M , also known as the rotation matrix, takes different
values for different crystal orientations and should be considered in separate
cases which are detailed in Appendix C. The simplest case, which we apply
to our nonlinear SPP structure is oriented along the (100) crystal orientation
41
with no rotation. In this simplifying case, together with the degenerate input
approximation, the nonlinear polarization (III.18) is simplified to
χdM :: E
aEbEc = χd(xˆE
a
xE
b
xE
c
x + yˆE
a
yE
b
yE
c
y + zˆE
a
zE
b
zE
c
z) (III.19)
In the following chapters we are going to use this simplest configuration to illustrate
nonlinear SPP properties without the cumbersome geometry distractions.
42
CHAPTER IV
HARMONIC GENERATION IN SURFACE PLASMON POLARITONS
Theoretical Development of Third Harmonic Generation of Nonlinear SPPs
This work was published in volume 696 of Frontiers in Optics (FiO)/Laser
Science XXVI (LS) conference in 2010. Dr. Miriam Deutsch set up the proposed
frame of the work I was the primary contributor to all the derivations and analysis
of the results.
In this chapter, we theoretically demonstrate that it is possible to utilize the
effect of third order susceptibility χ(3) in noble metals to induce third harmonic
generation (THG) from an SPP wave. The research on second harmonic generation
(SHG) on the metal surface by SPPs can be traced back more than two decades ago
on and it has been explored on a variety of metals, e.g. gold, silver and aluminum
[70] [71], as well as on different geometric configurations, e.g. bulk, grating an film
[71] [72] [73], all based on second order susceptibility χ(2). The effect of the intrinsic
third order susceptibility χ(3) of metals, however, were not covered until recently
by Novotny and et. al. [56] [67] in the SPP four wave mixing experiments.
We start our study under the non-depletion approximation as stated in
43
previous chapter. The wave equation (III.4) is
∇×∇× E+ ε
c2
∂2E
∂t2
= −µ0∂
2PNL
∂t2
(IV.1)
The steady-state solution of this wave equation has the form
E(r, t) = E(r)e−iωst (IV.2)
The spatial part of the field, E(r), can be written as the sum of the transverse
and longitudinal components
E = E⊥ + E‖ (IV.3)
where the transverse component relative to the propagation wavevector is labeled
by ⊥ and the longitudinal component by ‖ . Using the Helmholtz theorem in
Appendix B, the equation can be transformed according to:
∇×∇× E(r)⇒ −(ik)2E⊥(k) (IV.4)
while the temporal part, by Fourier transformation into frequency domain, is given
by
∂2
∂t2
E(r)⇒ −ω2sE(ω) (IV.5)
where ωs represents steady-state solution frequency. Now the nonlinear wave
equation(III.4) can be written as
−(ik)2E⊥(ω,k)− εω
2
s
c2
E(ω,k) = µ0ω
2
sP
NL(ωs,k) (IV.6)
44
The transverse ⊥ and longitudinal ‖ components of the equation can be
achieved by respectively applying k× and k· on (IV.6) which produces the two
separate equations:
−(ik)2E⊥ − εω
2
s
c2
E⊥ = µ0ω
2
sP
NL
⊥ (IV.7a)
−εε0k ·E‖ = k ·PNL‖ (IV.7b)
Performing an inverse Fourier transform on the equations above gives the
equations in component form [40]:
∇2E⊥(r, ω) + εω
2
s
c2
E⊥ = −µ0ω2sPNL⊥ (IV.8a)
∇ · [εε0E‖ +PNL‖ ] = 0 (IV.8b)
These are the central equations we are going to focus on.
To solve (IV.8a) and (IV.8b), first the explicit expression of PNL is needed.
Continuing from the discussion in last chapter, here a cubic medium (noble metal)
in (100) orientation with ϕ = 0 rotation is constructed for the geometry setup, so
that (III.18)
PNL = χ1122(E
b·Ec)Ea+χ1212(Ec·Ea)Eb+χ1221(Ea·Eb)Ec+χdM :: EaEbEc (IV.9)
where χd = χ
(3)
1111 − χ(3)1122 − χ(3)1212 − χ(3)1221, is now simplified to
Pi = 3χ
(3)
1212Ei(E ·E) + (χ(3)1111 − 3χ(3)1212)E3i (IV.10)
where i = x, y, z is the Cartesian component index.
45
We define two constant factors A ≡ χ(3)1111 and B ≡ 3χ(3)1212, and only include
terms of third harmonic generation:
Pi = BEi|E|2 + (A− B)E3i (IV.11)
The incident SPP wave form is predetermined as
E = (xˆEx0 + zˆEz0)e
−iω1t+ik1·r (IV.12)
The subscript 0 denotes amplitude value of the field at the interface, ω1 is the carrier
frequency and k1 = k(ω1) is the propagation wavevector of the SPP. Consequently
the components of PNL are written explicitly as
PNLx = BEx|E|2 + (A− B)E3x = Ex0(A|Ex0|2 +B|Ez0|2)e−iω3t+iq·r (IV.13a)
PNLy = 0 (IV.13b)
PNLz = BEz|E|2 + (A− B)E3z = Ez0(A|Ez0|2 +B|Ex0|2)e−iω3t+iq·r (IV.13c)
We the ω3 = 3ω1 frequency, the triple wavevector is defined q ≡ 3k1, as well
as the wavevector at triple frequency K ≡ k(ω3). Since the nonlinear polarization
has the dependency factor e−iω3t+iq·r, the steady-state frequency is determined by
ωs = ω3 = 3ω1, and the solution to (IV.8b) and (IV.8a) can be solved separately
in the dielectric and metal. Using the method of undetermined coefficients, and
considering the fact that PNL = 0 in the linear dielectric, the solutions take the
46
forms:
in dielectric: Ed = Ed0e
iKd·r (IV.14a)
in metal: Em = Em0e
iKm·r +Q0e
iq·r (IV.14b)
The subscript 0 denotes that the amplitude values are taken on the interface
between the dielectric d and the metal m. Here the Em0 term is a free wave
term that represents a radiating electro-magnetic field from the interface, while Q0
term is a driven term that represents an oscillating field acting as the secondary
source to the generated radiating wave. The interface amplitudes Ed0 and Em0, as
well as the wavevectors Kd, Km, which depend on the tripled frequency ω3, are
determined by boundary conditions. The interface amplitude Q0 = Q⊥ +Q‖, on
the other hand, is solved directly from the nonlinear polarization PNL. By solving
(IV.8a) and (IV.8b) we get
Q⊥ =
µ0ω
2
3P
NL
⊥
K2m − q2
(IV.15a)
Q‖ = −
PNL‖
εmε0
(IV.15b)
thus Q0 = Q⊥+Q‖ is determined. The longitudinal and transverse components of
a field are defined relative to the wavevector k according to Helmholtz theorem in
Appendix A. For example the polarization has the components defined as
P‖ =
P · k∗
|k|2 k (IV.16a)
P⊥ =
k∗ × (P× k)
|k|2 = P−P‖ (IV.16b)
47
For the next step, we need to determine all the wavevectors Kd, Km and
q. As stated above, q is already known and solely determined by form of the
nonlinear polarization PNL so that q = 3k1 where k1 = kspp(ω1) is the incident
SPP propagation wavevector at frequency ω1. Since we are looking for a steady-
state solution, all the wavevectors have to match along the x direction, so that
Kdx = Kmx = qx (IV.17)
With the relations that in the media K2dx+K
2
dz = εd(ω3)
ω23
c2
and K2mx+K
2
mz =
εm(ω3)
ω23
c2
, the z components of the wavevectors are now determined too. Thus
all the wavevectors in the 3ω frequency regime are determined at this point. The
imaginary part in Kd dominates, which indicates that the THG on the dielectric
side is still bound onto the surface, shown in Fig. 4.1. In the metal, the real part
in Km dominates, which indicates that the THG in the metal is able to radiate to
distance, shown in Fig. 4.2.
Next, to achieve all the field amplitudes, it is convenient to define the
wavevector angle θ by wavevector components:
tan θ ≡ kz
kx
(IV.18)
Also for the electric fields Ed0 and Em0, it is convenient to define an assistant
angle ϕ that ϕ = θ +
pi
2
because k ·E = 0 and it satisfies
tanϕ =
Ez
Ex
(IV.19)
48
Figure 4.1 A graph of the THG solution in vacuum. The THG solution in vacuum
is still bounded on the surface. Here the E-field component Ex is used in the graph.
Figure 4.2 A graph of the THG solution in metal. The THG solution in metal is
a propagating wave. The propagation angle is determined by the wavevector Kd.
Here the E-field component Ex is used in the graph.
Now we can find out all the field amplitudes on the interface by matching the
boundary conditions: The boundary condition for E along x direction at z = 0 is
xˆ · Ed = xˆ · Em + xˆ ·Q (IV.20)
49
from here we will talk about the field amplitudes on the interface and the subscript
0 is dropped for conciseness. This boundary condition leads to
Edx = Emx +Qx (IV.21)
or equivalently
Ed cosϕd = Em cosϕm +Q⊥ cosϕq +Q‖ cos θq (IV.22)
where the subscripts d, m and q denote the angles are calculated using wavevectors
Kd, Km and q respectively.
Similarly, the boundary condition for B along y direction at z = 0 is
yˆ ·Bd = yˆ ·Bm (IV.23)
which give
EdKd = EmKm + qQ⊥ (IV.24)
Together with the two equations (IV.22) and (IV.24) from boundary conditions,
the amplitude of electric fields on the interface can be solved
Ed =
Q⊥(Km cosϕq − q cosϕm) +KmQ‖ cos θq
Km cosϕa −Kd cosϕm (IV.25a)
Em =
Q⊥(Kd cosϕq − q cosϕd) +KdQ‖ cos θq
Km cosϕa −Kd cosϕm (IV.25b)
On the other hand, the Cartesian components of the electric fields on the interface
can be achieved using the boundary conditions while considering the conversion
50
that
Qx = Q‖ cos θq +Q⊥ cosϕq (IV.26a)
Qz = Q‖ sin θq +Q⊥ sinϕq (IV.26b)
then the equation (IV.24) is now
KdzEdx −KdxEdz = KmzEmx −KmxEmz + qzQx − qxQz (IV.27)
which leads to
Edx(tan θa − tanϕa) = Emx(tan θm − tanϕm) + (Qx tan θq −Qz) (IV.28)
The equations of the boundary conditions (IV.21) and (IV.28) can be solved to get
the x components of the fields:
Edx =
Qx(tan θq − tan θm + tanϕm)−Qz
(tan θd − tanϕd)− (tan θm − tanϕm) (IV.29a)
Emx =
Qx(tan θq − tan θd + tanϕd)−Qz
(tan θd − tanϕd)− (tan θm − tanϕm) (IV.29b)
And the z components of the electric field can be achieved via the relations that
Ez = Ex tanϕ accordingly. Up to this point, we achieved the solutions of the third
harmonic generation problem to the differential equations (IV.8a) and (IV.8b).
Conversion Efficiency and Propagation Angles
The propagation direction, the angle of propagation ϑ relative to the x axis,
of the free wave in the THG is determined by its wavevector Km via
tanϑ =
Re[Kmz]
Re[Kmx]
(IV.30)
51
When the imaginary part of Km is small, the two angles θ and ϑ are approximately
the same since |Km| ∼ Re[Km].
For quantitative illustration as well as for simplicity, we use the metal dielectric
function for silver given by the Drude model as an example,
ε(ω) = ε∞ −
ω2p
ω2 + iγpω
. (IV.31)
where for silver [54] the background dielectric constant ε∞ = 5.1, the plasma
frequency ωp = 8.65eV and the damping constant γp = 0.0045eV . Then the
expression of the propagation angle can be written explicitly in terms of the incident
frequency ω of the SPP wave:
tanϑ =
Kmz
Kmx
=
√
K2M − q2x
qx
=
√
K2M − (3kx(ω1))2
3kx(ω1)
=
√
ω23εm(ω3)− 9ω2
εm(ω1)
εd(ω1) + εm(ω1)
3ω
√
εm(ω1)
εd(ω1) + εm(ω1)
=
√√√√√√√ε∞ − ω2p9ω21 + 3iω1γp −
ε∞ −
ω2p
ω21 + iω1γp
εd(ω1) + ε∞ − ωp
ω21 + iω1γp√√√√√√√√
ε∞ −
ω2p
ω21 + iω1γp
εd(ω1) + ε∞ −
ω2p
ω21 + iω1γp
(IV.32)
here ω1 is the incident SPP wave frequency and ω3 = 3ω1 is the tripled frequency.
When the incident SPP frequency ranges within the optical range, e.g. 400nm ∼
700nm, the radiation angle of the THG ranges from 46◦ to 56◦, pointing off from
52
the interface into the metal, as indicated in Fig. 4.3.
Figure 4.3 A plot of the propagation angle. The propagation angle of the THG
v.s. the incident SPP wavelength. The angle is analytically determined by (IV.32).
With all the field strength relations available, we may now discuss the
conversion efficiency from the incident SPP wave to the propagating THG signal
wave. Here we choose the incident SPP to have a field amplitude of Ex = 10
4V/m,
which is proven to not violate the non-depletion approximation as we will see later.
This amplitude also corresponds to a 1W continuous laser focused onto a 1mm2
spot. The incident wavelength of the SPP is chosen to be 690nm for convenience
in later calculations.
First, we get the magnetic field for SPPs considering they are TM mode[74]
Hy =
By
µ0
= ε0c
2By =
ε0c
2
ω
(kzEx − kxEz), (IV.33)
The energy density of the SPPs are given in [74]. In the dielectric, there is no
53
dispersion, so on the surface of the interface,
ud =
1
4
[ε0E · E∗ + µ0H ·H∗] (IV.34)
where
E · E∗ = E0e−i(ωt−k·r) · E∗0e+i(ωt−k
∗·r) = (ExE
∗
x + EzE
∗
z )e
−2[Im(kx)x+Im(kz)z] (IV.35)
and
H ·H∗ = HyH∗ye−2[Im(kx)x+Im(kz)z] (IV.36)
The total energy (per unit area) associated with the surface polaritons is determined
by integration over z, the energy per unit surface area being
Ud =
∫ ∞
0
uddz =
ud
8Imkaz
=
1
8Imkaz
[ε0(ExE
∗
x + EzE
∗
z ) + µ0HyH
∗
y ] (IV.37)
For SPPs in a metal, the effective energy density on the surface is
um =
1
4
Re
[
d(ωε)
dω
]
ω
E(x) · E∗(x) + 1
4
µ0H(x) ·H∗(x) (IV.38)
And the total energy density in the metal (per unit area) is
Um =
∫ 0
−∞
umdz (IV.39)
The total energy over all space is
Utotal = Um + Ud (IV.40)
The third harmonic wave has the same expressions for the energy densities except
for a tripled frequency ω3.
54
Assuming the metal/dielectric interface is between silver and vacuum such
that εd = 1, then the conversion efficiency from SPP to THG wave is
ηeff =
Um(ω3)
Utotal(ω1)
∼ 0.1% (IV.41)
This result is consistent with the non-depletion approximation that the energy
drained from SPP wave by THG signal is small enough that the SPP wave remains
relatively unaltered.
Adjustment to Realistic Noble Metals
According to the Drude model, silver is absorptive at wavelength λ = 690nm,
but it would be transparent when frequency is tripled at wavelength λ = 230nm,
which is known as the ultraviolet transparency for metals. However, the noble
metals have increasingly strong absorption due to interband transitions, as seen in
Fig. 2.1b.
For the case of silver, however, we found in literature data a very narrow
transparency window just above the plasma frequency in the wavelength range of
293nm ∼ 314nm, so that within this wavelength, the real part of the dielectric
function is positive (transparent) while the imaginary part (absorption) is not
dominating so a portion of the generated THG wave can still penetrate and
propagate [37]. In order to have the THG wave result in the transparency
window, the incident SPP wavelength should be tuned to the infrared range
55
of 880nm ∼ 941nm. This corresponding frequency range of the SPPs and the
resultant transparency window frequency range are shown in Fig. 4.4.
Figure 4.4 Drude model vs realistic metals. Close agreement is shown between
Drude model and real data for the dielectric function of silver. The shaded
green represents the transparency window and the (dense) shaded pink represent
corresponding suggested wavelength of the incident SPPs. The literature data is
cited from [37].
To reduce the absorption of the THG wave in silver, a metal film can be used.
When a silver film of 50nm thickness is used, the generated THG wave has 35% of
the amplitude (10% of the intensity) transmitted through the metal film.
The interband transition of silver also alters the calculation in the propagation
angle of the THG in metal, via the realistic data of εm(ω). The ϑ− λspp curve has
an obvious red shift but keeps the basic shape similar to the analytical one in Drude
56
model. The resultant THG wave can be easily separated from the incident SPPs
since the radiation direction of the THG wave clearly steers away at a rather large
angle relative to the SPPs, as shown in the shaded area in Fig. 4.5.
Figure 4.5 A plot of propagation angle with silver data. The angle of propagation
for the THG in metal is altered from the case with the Drude model. The solid
line is the angle predicted with the Drude model and the dotted curve is angles
calculated from literature data of silver. The shaded area indicates the suggested
wavelength for the incident SPP. Using the suggested wavelength, the resultant
THG will fall into the transparency window of silver and the radiation angle is
shown in dots in the graph.
Conclusions
This chapter demonstrated that a third harmonic generation (THG) wave is
induced by SPPs in a nonlinear noble metal. The energy conversion efficiency is
57
calculated to be about 0.1% with the intensity of a common used laser in labs,
e.g. with amplitude E ∼ 104V/m. When the frequency of the induced THG wave
is above the plasma frequency of the metal, the THG wave is predicted to radiate
through the metal if the Drude model is adapted for optical properties of the metal.
For realistic metals, however, the radiation of the THG wave in the metal is strongly
absorbed at and above the plasma frequency. A special case among the noble metals
is silver. Silver, according to the experimental data in the literature, has a narrow
ultraviolet transparency window around the plasma frequency. Our calculation for
silver shows reasonable transmittance of the THG wave through a 50nm film and
the direction of THG radiation is clearly separated from the incident SPPs by its
propagation angle.
58
CHAPTER V
PLASMONIC SOLITONS
Introduction to Optical Solitons
This work was unpublished. Dr. Miriam Deutsch set up the proposed frame
of the work I was the primary contributor to all the derivations and analysis of the
results.
A solitary wave, often referred to as a soliton, describes a special kind of
nonlinear wave, which keeps its wave properties unchanged, such as the pulse shape
or beam diameter, during propagation in a medium. This kind of stable solitary
wave can exist only with the presence of nonlinearities in the host medium. The
spatial properties of a wave are changed because of the wave nature: a wave always
diffracts during propagation. The temporal properties of a wave are changed, such
as pulse broadening, if there is any dispersion present in the host dispersive medium.
The spatial diffraction and temporal broadening, however, can be compensated in a
medium with nonlinear response. Thus it is possible to achieve a stable wave form,
a soliton, when the broadening and the diffraction are completely counteracted by
nonlinear effects [75]. The fundamental mathematical models of solitons were built
59
in the 1960s [76] [77], and the concept of solitons has since evolved to a much
broader range in different areas of physics [78][79][80][81][82][83].
In the context of nonlinear optics, one very important example of optical
solitons, which is widely used in the field of modern communication technology, is
a stable optical pulse guided in an optical fiber [84]. This stable optical pulse is
formed by the balance of the (nonlinear) optical Kerr effect in the fiber and the
material’s (linear) chromatic dispersion [85][86]. The optical Kerr effect is a direct
result from the third order susceptibility χ(3) of the glass, and this optical Kerr
effect in the fiber is represented by the change of the optical refractive index ∆n,
which is proportional to the light intensity, I.
Optical solitons are classified as being temporal or spatial according to the
confinement types of the optical waves. For light with temporal confinement such
as an optical pulse, the chromatic dispersion, e.g. high frequency wave components
travel faster than low frequency components, introduces a frequency chirp and
therefore leads to broadening in the pulse. Due to different phase velocities of
frequencies in the optical pulse, frequency components are redistributed such that
higher frequencies are in the leading half of the pulse and lower frequencies in the
trailing half. In a nonlinear medium, this chirp may be compensated by the optical
Kerr effect, because the nonlinear interaction between the electric field and the
medium causes a down-shift of the frequencies in the leading half of the pulse, and
an up-shift in the trailing half. These two opposing effects, in some special cases,
60
can be balanced and thus the pulse keeps its amplitude envelope unchanged in time
and this is called an optical temporal soliton, as shown in Fig. 5.1.
For an optical beam in space, an analogy to the pulse can be used, where
the confinement now is in one or more transverse spatial dimensions instead of
in time for an optical pulse. A light beam is usually diffracted, spreading the
beam diameter during the wave propagation. The nonlinear refractive index here
acts as a focusing lens and induces a self-focusing (or self-defocusing) effect with a
positive (or negative) Kerr coefficient [83]. When the self-focusing counteracts and
cancels the effect of diffraction, the light beam may maintain its diameter during
propagation in space, and it is called an optical spatial soliton as shown in Fig. 5.2
[40].
(a) (b)
Figure 5.1 A schematic of a temporal soliton formed by the balance between linear
dispersion and optical Kerr effect. (a) The optical pulse is chirped and broadened
by dispersion in a linear medium. (b) The dispersion can be balanced to form an
unchanged wave packet by the optical Kerr effect.
61
Figure 5.2 A schematic of a spatial soliton formed by the balance between the
diffraction and self-focusing. (a) The optical beam diffracts in a linear medium
while propagating along x−axis. (b) The diffraction is counteracted by the self-
focusing effect since the refractive index is higher towards the center of the optical
beam. (Figure taken from ref. [83])
The first experimental observation of optical temporal solitons were discovered
in optical fibers in 1980 by Mollenauer et al. using pico-second laser pulses [84].
When an optical pulse propagates in an optical fiber, the chromatic dispersion
induces a group-velocity dispersion (GVD) and the nonlinear refractive index leads
to a self-phase modulation (SPM). The GVD is the group velocity dependency on
angular frequencies, namely the derivative of the inverse group velocity, vg, with
respect to the angular frequency, ω:
GVD ≡ ∂
∂ω
1
vg
=
∂2β
∂ω2
(V.1)
where β is the propagation wavevector. The SPM is a nonlinear phase delay, ∆φ(t),
which arises from the varying optical refractive index ∆n(I), where I is the optical
62
intensity. In an optical pulse, this phase delay leads to a change of the instantaneous
frequency, ω(t), according to ∆ω(t) =
d∆φ(t)
dt
. Consequently the change of the
instantaneous frequency is relative to the pulse instantaneous intensity by [87]:
∆ω(t) ∼ −∂∆n(I(r, t))
∂t
(V.2)
For a medium with optical Kerr effect such that ∆n ∝ I, (V.2) gives ∆ω(t) < 0
for the leading half of the pulse where the instantaneous intensity I(t) is ramping
up, and ∆ω(t) > 0 for the trailing half. At the communication wavelength, 1.55µm,
the silica-glass fiber has an anomalous GVD (GVD< 0) and a positive Kerr effect
(χ(3) > 0). An anomalous GVD means the group velocity is faster for higher
frequencies and slower for lower frequencies in an optical pulse. The anomalous
GVD gives the optical pulse a down-chirp, namely the instantaneous frequency
of the pulse decreases with time, as indicated in Fig. 5.1(a). At the same time,
according to (V.2), the SPM will tune down frequency when the intensity increases,
which is the case for the leading half of the optical pulse, and tune up frequency
when the intensity decreases, which is the case for the the trailing half. With a
carefully chosen pulse envelope (a sech function), the two competing effects over
the frequencies can completely cancel each other and a temporal soliton is achieved,
as shown in Fig. 5.1(b).
A spatial soliton is formed in a scenario similar to a temporal soliton. On
one hand, the wave nature of an optical beam tends to diffract the beam diameter
during propagation in a medium, as indicated in Fig. 5.2(a). On the other hand,
63
with the presence of the positive Kerr effect (χ(3) > 0), the center of the beam,
which of high optical intensity, a induces larger optical refractive index than the
rim area of the beam. This index inhomogeneity works in a similar way to a convex
lens that focuses a light beam, thus results a self-focusing effect of the light beam,
as shown in Fig. 5.2(b). For a beam with a properly chosen initial shape and
intensity, the self-focusing of the beam may exactly counteract the diffraction so
that the optical pulse remains a confined, self-trapped beam. Thus a spatial soliton
is achieved [88] [89].
The traditional wave-packet (temporal) solitons discussed in previous sections
are called bright solitons since they are positive energy packages against background
of zero field intensity, as shown in Fig. 5.3a. When the GVD is normal (GVD> 0),
however, a bright soliton is not supported by the nonlinear optical fibers discussed
above. Instead, a special “dark” pulse, which exists as an absence of energy on a
continuous wave (CW) background, turns out to be stable and thus called a dark
soliton, as shown in Fig. 5.3b. With a positive Kerr effect (χ(3) > 0), the bright
temporal solitons can be achieved in a nonlinear medium with an anomalous GVD
(GVD< 0) while that the dark temporal solitons can be achieved with a normal
GVD (GVD> 0). For the case of spatial optical beams, diffraction plays a role
similar to the GVD in the temporal soliton case: the self-focusing nonlinearity from
positive Kerr effect (χ(3) > 0) gives bright spatial solitons and the self-defocusing
nonlinearity from negative Kerr effect (χ(3) < 0) gives dark spatial solitons.
64
(a) Bright soliton (b) Dark soliton
Figure 5.3 A bright soliton and a dark soliton. (a) A bright temporal soliton is
an optical wave packet which maintains the envelope against zero background. It
can be achieved when the nonlinear medium has an anomalous GVD (GVD< 0).
The envelope function for the intensity envelope function of the bright soliton is a
|sech|2 function. (b) A dark temporal soliton is a stable “dark” hole which keeps its
shape against non-zero background. It can be achieved when the nonlinear medium
has a normal GVD. The intensity of a dark soliton is a |tanh|2 function.
The general theory on solitons has been developed long before the observation
of the optical solitons and applied to a variety of different nonlinear systems [76]
[90] [77]. There are a number of partial differential equations (PDEs) with soliton
solutions in different forms and a few may be solved exactly, such as the Kortveg-de
Vries (KdV) equation, the sine-Gordon equation and the nonlinear Schro¨dinger
(NLS) equation [75][91]. For a light-wave packet traveling in a waveguide, the
evolution of the envelope is commonly described by a NLS equation [86], which
65
results in an optical temporal soliton solution. The general solutions of the NLS
equation can be achieved by the inverse scattering transform [92], a method for
nonlinear systems similar to the Fourier transform in linear systems and the
solutions are known to be solitons.
In this dissertation, we discuss nonlinear SPPs guided by a planar structure,
which is an interface between a semi-infinite metal with Kerr nonlinearities and a
semi-infinite linear dielectric. We focus on the temporal equation for the envelope
of an SPP pulse, and find the conditions for a soliton-type solution. Early works
on spatial SPP solitons are discussed in detail by pioneer researchers in [83], and
the solitons on the surface comprised of a nonlinear dielectric and a linear metal
are discussed in [93]. The early works about SPP solitons focused on the dielectric
nonlinearity but by and large ignores the metal nonlinearity. However, as we stated
in previous chapters, the metal nonlinearity is comparable to or even greater than
the nonlinearity in many dielectrics such as glass. With the negative permittivity
of a metal, nonlinearities in the metal may lead to different phenomena than
nonlinearities in a dielectric. In this chapter, our research is presented to reveal the
soliton properties while the metal nonlinearity is taken into account.
The Temporal Nonlinear Schro¨dinger Equation of SPPs
The linear SPP is a two-dimensional guided wave propagating at the interface
66
between a dielectric and a metal. In this section we study the new features of the
nonlinear SPPs when the metal nonlinearity is taken into account.
We consider only the third order nonlinear susceptibility χ(3) in the metal,
but ignore its sum frequency effect, such as was discussed in Ch.IV in the sections
concerning THG, because the sum frequency effect occurs at very low efficiency.
Therefore the metal is regarded as a Kerr medium and the nonlinear polarization
PNL in (III.5) has the same frequency as the original incident light E
PNL(ω) = ε0χ
(3)(ω = ω + ω − ω)E(ω)E(ω)E(ω) (V.3)
Combining the time derivatives in the nonlinear wave equation(III.4) gives
∇×∇×E+ 1
c2
∂2
∂t2
εME = 0 (V.4)
where εM = εm+εNL is the total dielectric function of the metal, including both the
linear contribution εm = 1+χ
(1), and the nonlinear contribution εNL = χ
(3)|E|2, at
response frequency ω. Here we prefer using the optical dielectric function ε of the
metal instead of the optical refractive index n even when n is commonly used in the
analysis for nonlinear fiber optics [86] [94]. We do this because ε can be negative
or complex in our discussion and the traditional definition of n =
√
ε often leads
to inconvenience and confusion.
The divergence free equation ∇·E = 0 holds inside the isotropic metal except
at the interface, such that ∇ × ∇ × E = −∇2E. Consequently the Cartesian
components are not coupled in (V.4) with the chosen coordinates shown in Fig. 2.2
67
and (V.4) can be written in component form:
xˆ : − (∂2x + ∂2y + ∂2z )Ex +
1
c2
∂2
∂t2
εMEx = 0 (V.5a)
zˆ : − (∂2x + ∂2y + ∂2z )Ez +
1
c2
∂2
∂t2
εMEz = 0 (V.5b)
where ε0µ0 = 1/c
2. Considering the TM nature of linear SPP profile, Ex and Ez
have a fixed relation Ex = ζEz where ζ is a piecewise constant in the dielectric
and in the metal. Since the two equations are linearly dependent, we may consider
Ez first under the scalar approximation to eliminate the complexity of the vector
nature of the nonlinear wave equation. The equation derivation process for Ex is
identical to Ez. Apply Fourier transform
Ez(r, t) =
1√
2pi
∫
E˜z(r, ω − ω0)e−i(ω−ω0)tdt (V.6)
upon equation (V.5b), with ω0 being the central frequency of the incident SPP
wave, and it gives
∇2E˜z + ω
2
c2
εME˜z = 0 (V.7)
The tilde “∼” represents that the quantity is a function of frequency ω instead
of time t. Separating the propagation factor eiβ0x from E˜z (where β0 being the
propagation constant to be determined from the eigensolution of the equation) the
electric field E˜z can be written as
E˜z(x, y, z, ω − ω0) = Fm(y, z)A˜z(x, ω − ω0)eiβ0x (V.8)
Here Fm(y, z) and A˜z(x, ω − ω0) are the transverse and longitudinal functions
68
respectively. The differential equation (V.7) is now
A˜z(∂
2
y + ∂
2
z )Fm + Fm(∂
2
xA˜z + 2iβ0∂xA˜z − β20A˜z) +
ω2
c2
εMFmA˜z = 0 (V.9)
For a long, smooth optical pulse the slowly varying amplitude approximation
|∂2xA˜z| ≪ |β0∂xA˜z| is satisfied so that the higher order term ∂2xA˜z ≃ 0 is ignored.
Applying the standard method of variable separation to the partial differential
equation (V.9) yields
(∂2y + ∂
2
z )Fm
Fm
+ εM
ω2
c2
= −2iβ0∂xA˜z − β
2
0A˜z
A˜z
≡ β˜2 (V.10)
where β˜ is an undetermined coefficient and it can be a function of ω.
The differential equation for Ez in the dielectric medium has a similar
expression:
(∂2y + ∂
2
z )Fd
Fd
+ εd
ω2
c2
= −2iβ0∂xA˜z − β
2
0A˜z
A˜z
≡ β˜2 (V.11)
with A˜z, β0 and β˜ the same as required by the boundary conditions, but with a
different transverse functions Fd and a different dielectric function εd. We can
make the further assumptions that the SPP wave extends to ±∞ along the y
axis, and that the derivatives with respect to y are negligible Therefore we have a
reduced 2-dimensional problem with x and z dependencies only. Consequently Fd
and Fm are functions of z only and the partial derivatives against z become total
derivatives. The two equations (V.10) and (V.11) on different sides of the interface
69
can be combined by writing F and εT as piecewise functions
F (z) ≡

Fd(z) z > 0 (Dielectric)
Fm(z) z < 0 (Metal)
(V.12)
εT (z, ω) ≡

εd(ω) z > 0 (Dielectric)
εM(ω) z < 0 (Metal)
(V.13)
Here εM is the nonlinear dielectric function of the metal and possesses the property
that ReεM < 0 and |ReεM | ≫ |ImεM |. The parameter εd is the dielectric
constant, which is close to unity and |εM | ≫ |εd| & 1. For partial differential
equations, the standard procedure to achieve the common coefficient β˜ is to solve
the eigenfunctions:
d2F
dz2
−
[
β˜2 − εT ω
2
c2
]
F = 0 (V.14a)
2iβ0
∂A˜z
∂x
+ (β˜2 − β20)A˜z = 0 (V.14b)
We obtain the linear solution of function F by letting εM = εm and use it as the
first order approximation for the nonlinear case to achieve the function A˜z and thus
the first equation in (V.14) is written explicitly as a piecewise eigenequation
d2F
dz2
=

(
β˜2 − εdω
2
c2
)
F z ≥ 0(
β˜2 − εmω
2
c2
)
F z < 0
(V.15)
where εm is the linear dielectric function for the metal.
70
Now consider the physical restrictions that the waves from both sides match
at the boundary and decay at ∞. When εm < 0, this is possible only if β˜2 > εdω
2
c2
.
Thus the boundary condition for Ez gives the transverse solution:
F (z) =

εme
−z
√
β˜2−εd
ω2
c2 z ≥ 0
εde
+z
√
β˜2−εm
ω2
c2 z < 0
(V.16)
which is exactly the linear solution for SPP stated in previous chapters.
The boundary condition for Ex, on the other hand, leads to
dF
dz
∣∣∣∣
z=0+
=
dF
dz
∣∣∣∣
z=0−
(V.17)
resulting again in the linear SPP dispersion relation
β˜2 =
ω2
c2
εmεd
εm + εd
≡ β˜2L (V.18)
For the the nonlinear case, an ansatz is achieved by perturbing β˜, to the first
order, from the linear β˜L such that now β˜ ≡ β˜L +∆β˜. The two photon absorption
effect of bulk noble metals is not fully studied but known to be very weak[95], and
it is not directly related to our topic. To keep the physics simple and succinct,
we ignore the two photon absorptions related to the χ(3) and keep χ(3) a pure real
number. The discussion, however, can be extended to a χ(3) with both real and
imaginary parts with nontrivial conversions [96].
It is useful to Taylor expand β˜L about the carrier frequency ω0
β˜L(ω) = β0 + (ω − ω0)β1 + 1
2
(ω − ω0)2β2 + 1
6
(ω − ω0)3β3 + · · · (V.19)
71
and keep terms up to the second order, where the constant β0 = β˜L|ω=ω0 is the
propagation constant at carrier frequency ω0 and the other coefficient constants are
defined as
βn ≡ d
nβ˜L
dωn
∣∣∣∣∣
ω=ω0
(V.20)
where β1 = 1/vg is the inverse of the group velocity, and β2 is known as the group
velocity dispersion, or GVD.
With approximation to the first order, the right hand side of (V.14b) is
β˜2 − β20 ⋍ 2β0(β˜ − β0) (V.21)
and the equation is written as
∂A˜z
∂x
= i
[
β˜L(ω) + ∆β˜ − β0
]
A˜z (V.22)
Now using the Taylor expansion of β˜L and keeping second order accuracy,
∂A˜z
∂x
= i
[
(ω − ω0)β1 + 1
2
(ω − ω0)2β2 +∆β˜
]
A˜z (V.23)
Using an inverse Fourier transform
Az(x, t) =
1√
2pi
∫ +∞
−∞
A˜z(x, ω − ω0)e−i(ω−ω0)tdω (V.24)
on the equation above leads to the nonlinear schro¨dinger equation (NLS) equation
for the longitudinal slowly varying amplitude A
∂Az
∂x
= −β1∂Az
∂t
− i
2
β2
∂2Az
∂t2
+ i∆βAz (V.25)
72
A similar equation can be derived for the Ax component, too. Using the linear
relation that Ex = ζEz, the combination of the equations for Az and Ax gives the
NLS equation for the amplitude A
∂A
∂x
= −β1∂A
∂t
− i
2
β2
∂2A
∂t2
+ i∆βA (V.26)
where A = |A| and A = zˆAz + xˆAx.
This is the equation that governs the evolution of the amplitude of the SPP,
taking into account the Kerr effect of the supporting metal. It is the starting point
of plasmonic solitons and gives temporal soliton solutions. The coefficients β1 and
β2 are approximated to the first order by derivatives of the linear β˜L with respect
to angular frequency. ∆β is the perturbation to the linear wave factor with the
presence of χ(3) nonlinearity. It is under further consideration in next section. In
addition the attributes of the soliton solutions are discussed in detail in the following
sections.
Nonlinear SPP Dispersion Relations
In this section, we are going to determine the general properties of the
coefficients β1, β2 and ∆β in the NLS equation (V.26) from the last section. For
a TM wave guided along the x direction with propagation factor e−i(ωt−β˜x), the
nonlinear dispersion relation is derived from the field relation with the nonlinear
73
dielectric function, εM , in the nonlinear metal [83] [97] [98]:
dEx
dz
= − i
β˜
(
ω2
c2
εM − β˜2)Emz (V.27a)
d
dz
(εMEmz) = −iβ˜εMEx (V.27b)
Hy = −ωε0εM
β˜
Emz (V.27c)
We can get a first integral from the above relation equations:
(
dEx
dz
)2
=
ω2
c2
[
(εM − c
2β˜
ω2
)E2mz − εME2x − χ(3)(ExEmz)2 +
χ(3)
2
(E4x + E
4
mz)
]
(V.28)
In the next step we aim to obtain the nonlinear SPP dispersion relation. Using
(V.27a) for the left-hand side in equation (V.28) above:
(1− εMc
2
β˜2ω2
)E2mz = (εM−
c2β˜
ω2
)E2mz−εME2x−χ(3)(ExEmz)2+
χ(3)
2
(E4x+E
4
mz) (V.29)
To eliminate εM , the transverse profile and boundary conditions of the E field
are applied. In the dielectric, the transverse dependency is E(z) ∼ e−κdz where
the decay factor κd =
√
β˜2 − εdω
2
c2
, while the subscript d labels components in the
dielectric. Then the divergence free equation for the displacement field ∇ ·D = 0
gives Edz = iβ˜Ex/κd.
Considering components of the electric displacement, Ddz = Dmz, leads to
εdEdz = εMEmz . On the interface it gives
εM = iεd
β˜
κd
Ex
Emz
(V.30)
74
Using the expression for εM , (V.29) is now
χ(3)
2
E4x −
(
εM +
cεd
ωκd
)
E2x − i
εdβ˜
κ
ExEmz − χ
(3)
2
E4mz = 0 (V.31)
where all the E fields here take the values at the surface (z = 0). This equation
determines the nonlinear SPP dispersion relation since β˜ is implicitly determined
by ω,Ex and Emz . It is easy and straightforward to verify that it complies with
the linear dispersion relation by letting χ(3) = 0 and εM = εm.
For a general nonlinear dispersive metal, an additional relation may be required
between Ex and Emz . We first consider an ideal case that the losses in the metal
are small and negligible. This gives a pi/2 phase shift between Emz and Ex, thus
we can define a real-valued total amplitude ET > 0 with
E2T = |ET |2 ≡ |Emz|2 + |Ex|2 = E2mz − E2x (V.32)
since the total filed has the form zˆEmz + ixˆEx. Together with
εM = i
εdβ˜
κd
Ex
Emz
(V.33)
the first two relations are achieved by squaring (V.33), and gives the third relation
is given by multiplying (V.33) with E2mz :
E2x = −E2T
(εM)
2κ2d
ε2dβ˜
2 + (εM)2κ2d
(V.34a)
E2mz = E
2
T
ε2dβ˜
2
ε2dβ˜
2 + (εM)2κ2d
(V.34b)
ExEmz = −iE2T
εdεM β˜κd
ε2dβ˜
2 + (εM)2κ
2
d
(V.34c)
75
Finally we express all field components in terms of E2T . We use κ
2
d = β˜
2 − εdω
2
c2
in
(V.31) to solve for the nonlinear dispersion relation of β˜:
β˜2 =
ω2
c2
εMεd
(
εM − εd − E
2
Tχ
(3)
2
)
[(εM)2 − ε2d]−
E2Tχ
(3)
2
[(εM)2 + ε2d]
(V.35)
The nonlinear dielectric function of the metal εM = εm+χ
(3)E2T can be inserted
into (V.35) for a more explicit expression. Care must be taken when one tries to
achieve β˜ and operates a square root on the right hand side, because the square
root operation requires a proper cut in the complex plain so that the condition
Reβ˜ > 0 holds [99]. Using a Taylor expansion identity with respect to the small
quantity δ = χ(3)E2T :
(εm + δ)
√
a− δ
gx+ b
=
√
a
b
[
εm +
(
1− εm
2a
− gεm
2b
)
δ
]
+O(δ2) (V.36)
where the intermediate notations a = 2(εd−εm), b = 2εm(ε2m+ε2d) and g = 3ε2d−5ε2m
are used. To the first order of χ(3)E2T we have
β˜ =
ω
c
√
εdεm
εm + εd
− ω
c
χ(3)E2T
εd
4(εm + εd)
√
εd
εm(εm + εd)
(V.37)
The expression of εm is often known from tabulated data, or an analytic model,
e.g. Drude model as discussed in Ch.II, and we can calculate
∆β ≡ (β˜ − β˜L)|ω=ω0 =
ω0
c
χ(3)E2T
√
ε3m
16εm(εm + εd)3
∣∣∣∣∣
ω=ω0
(V.38)
For simplicity, we define the nonlinear factor
ν =
ω0
c
χ(3)
√
ε3m
16εm(εm + εd)3
∣∣∣∣∣
ω=ω0
(V.39)
76
Considering that ET ∼ Ae−iωteiβx, the NLS equation (V.26) is now written as
∂A
∂x
= −β1∂A
∂t
− i
2
β2
∂2A
∂t2
+ iν|A|2A (V.40)
The other coefficients β1 and β2 are given by evaluating to the first order using the
linear dispersion relation at ω = ω0
β1 ≡ dβ˜L
dω
∣∣∣∣∣
ω=ω0
and β2 ≡ d
2β˜L
dω2
∣∣∣∣∣
ω=ω0
(V.41)
Analytical and Numerical Solutions of Dark Solitons
Dark Soliton Solutions
In a lossless metal, the expression for εm from the Drude model is
εm =
ω2p
ω2
− ε∞ (V.42)
where ε∞ is the background dielectric constant at low frequency. All coefficients in
the NLS equation are real quantities and explicitly they are written as functions of
frequency ω:
β˜L =
ω
c
√
εd
√
ω2p/ω
2 − ε∞
ω2p/ω
2 − ε∞ − εd (V.43a)
∆βL =
ω0
c
χ(3)E2T
√
ε3d
16(ω2p/ω
2 − ε∞)(ω2p/ω2 − ε∞ − εd)3
(V.43b)
It is clear that ∆βL > 0 if χ
(3) > 0, which is consistent with the fact that the
larger (nonlinear) refractive index leads to a bigger wavevector. The β1 and β2
are calculated and graphed as in Fig. 5.4. Note that below the plasma frequency
77
ωp, the quantities β1, β2 and ∆β are all positive real numbers with lossless Drude
model (V.42).
ΩΒ1
Ω2 Β2
ΒL
Ω/c
0.0 0.2 0.4 0.6 0.8
ِΩSPP0
1
2
3
4
kSPP
Figure 5.4 The coefficients in the NLS equation. The coefficients β1 and β2 in the
NLS equation are calculated from the linear SPP dispersion relation using lossless
Drude model. β1 and β2 are positive below plasma frequency ωp. The vertical
dashed line is the asymptote of β˜L, which represents the frequency limit of a linear
SPP wave.
To solve the NLS equation (V.40) analytically, we transform A(x, t) =
78
A(ξ(x), τ(x, t)) by defining
τ = (t− β1x)/T0 (V.44a)
ξ = x/LD (V.44b)
u = A
√
|ν|LD (V.44c)
where T0 is the duration of the incident SPP pulse, LD = T
2
0 /|β2| is the SPP
dispersion length due to GVD and |ν| is the absolute value of ν. Then we have a
NLS equation with standard form
i
∂u
∂ξ
− 1
2
β2
|β2|
∂2u
∂τ 2
+
ν
|ν| |u|
2u = 0 (V.45)
For SPPs on noble metal surface, χ(3) ≡ Re[χ(3)] > 0 so that ν > 0. Also, with
normal group-velocity dispersion GVD≡ ∂
2β
∂ω2
> 0 so that → β2 > 0. We have now
i
∂u
∂ξ
− 1
2
∂2u
∂τ 2
+ |u|2u = 0 (V.46)
and this standard form of a NLS equation has a dark temporal soliton solution [86].
It is conventional to use a slightly different standard form with a positive second
order derivative term [92] [86]. By defining X = −ξ and U(X, τ) = u(ξ, τ), the
NLS equation turns to
i
∂U
∂X
+
1
2
∂2U
∂τ 2
− |U |2U = 0 (V.47)
The most general form of the solution is
U(X, τ) = U0B2 tanh[U0B2(τ − B1U0X)] + iB1e−iV 20 X+iθ0 (V.48)
79
where B1, B2 and θ0 are arbitrary constants which depend on initial conditions. For
a dark soliton, the boundary condition is
|U |ξ→±∞ = U0 (V.49)
where U0 is the background wave amplitude.
The coefficients B1 and B2 have a constraint that B
2
1 +B2 = 1. Thus we may
write B1 = sinφ and B2 = cosφ where the variable φ represents the “darkness”
of the soliton. For example when sin φ = 0 the minimum intensity of the solution
goes to 0 and we call it a “dark” soliton. When sin φ 6= 0, the minimum intensity
of the soliton is greater than 0 and the solution forms a “gray” soliton.
Switching back to the lab frame of reference with parameters t and x , the
solution is
A(x, t) =
A0
{
cosφ tanh
[
A0
√
ν cos φ√
β2
t−
(√
νβ1 −A0ν sinφ√
β2
)
A0 cosφx
]
+ i sin φ
}
eiA
2
0νx+iθ0
(V.50)
And the normalized intensity is
∣∣∣∣A(x, t)A0
∣∣∣∣2 = 1− cos2 φ sech2 [A0√ν cosφ√β2 t− A0 cosφ
(√
νβ1√
β2
+ A0ν sinφ
)
x
]
(V.51)
To satisfy this eigenfunction, the pulse duration, pulse length and group velocity
80
are required to be
τ0 =
√
β2
A0
√
ν cosφ
(V.52a)
L0 =
√
β2
(
√
νβ1 − A0ν sin φ
√
β2)A0 cosφ
(V.52b)
vg =
1
β1 −A0
√
ν sinφ
√
β2
(V.52c)
To illustrate the core physics, the simplest case is chosen, i.e. φ = 0 and θ0 = 0 in
(V.50) so the dark soliton has an amplitude and intensity equations
A(x, t) = A0 tanh
[
A0
√
ν√
β2
(t− β1x)
]
eiA0νx (V.53a)
|A|2 = |A0|2
{
1− sech2
[
A0
√
ν√
β2
(t− β1x)
]}
(V.53b)
Here A0 is the peak amplitude of the optical soliton. For the soliton, the peak
amplitude is not a free parameter but determined by the duration τ0 and other
constants as in (V.52a). In fact, (V.52a) determines the fundamental amplitude of
a soliton for a given duration τ0 by
A0 =
√
β2
τ0
√
ν cosφ
. (V.54)
Given a pulse duration, if the initial peak amplitude is slightly different to
the fundamental amplitude in (V.54), then the pulse will adjust itself during
propagation and eventually reach the stable soliton amplitude given by (V.54).
If the initial peak amplitude is very low, the nonlinear effect can not counteract
the linear dispersion and the pulse will be simply linearly dispersed. If the initial
81
peak amplitude is very high, complicated behaviors such as high-order soliton and
soliton break-up could occur [92] [91] [75].
Numerical Analysis with and without Metal Loss
The existence of the dark SPP soliton for the nonlinear SPP is demonstrated
by the numerical solution to the NLS equation (V.46) where the amplitude of the
electric field is already normalized to the stable soliton solution amplitude A0, and
in the initial condition we also set the arbitrary phase φ = 0. As shown in Fig. 5.5,
where the loss in the metal is ignored, the SPP pulse propagates along the x-
direction on the metal surface while its temporal intensity profile is conserved over
distance.
When the amplitude-duration relation in (V.54) is not satisfied, the optical
pulse is not an eigensolution to the NLS equation. We show numerically it will
either break down to separate pulses for an initial pulse with larger amplitude,
as shown in Fig. 5.6a, or oscillate during propagation for a pulse with smaller
amplitude, as shown in Fig. 5.6b. In the case of larger amplitude, two small “gray”
pulses appear symmetrically on each side of the main “dark” pulse, with the main
pulse width narrowed down. In the case of small amplitude, small continuous
oscillations propagate out of the center pulse on both sides. The initial shape of
the fundamental “dark” pulse is a tanh function, which is an odd function, or
odd symmetry. Therefore, the main pulse keeps “dark”, i.e., center intensity being
82
Figure 5.5 The stability of the soliton presented by numerical calculations. The
temporal profile of the SPP pulse does not change while propagating along the
x-direction on the nonlinear metal surface. Thus a stable dark soliton is achieved.
zero, regardless of the initial conditions. For reference numerical values, the SPP
related material constants used in the numerical calculation are listed in Tab.D.1
(for vacuum/silver) [41].
In our numerical solution for the lossless metal model, an SPP pulse with 10ps
duration is chosen at the communication wavelength 1.55µm. The peak amplitude
is calculated by the threshold equation (V.52a), which may be attained using a
83
(a) (b)
Figure 5.6 A schematic of pulse propagation when the soliton condition is not
satisfied. (a) The optical pulse breaks up when the initial amplitude is greater
than the threshold. [100] (b) The optical pulse has an oscillation when the initial
amplitude is less than the threshold.
pulsed laser at an average power of 0.5W with 80 MHz repetition value. The pulse
evolution for different initial amplitudes is shown in Fig. 5.5 and Fig. 5.6, and the
numerical parameters are shown in Tab.D.1.
We now turn to the more realistic case of a noble metal with Ohmic loss. The
noble metal has finite absorption in optical range, thus a linear lossy term is added
into the lossless NLS equation:
i
∂U
∂X
+
1
2
∂2U
∂τ 2
− |U |2U = iΓU (V.55)
where Γ ≡ LD/(2β1) is a linear loss parameter representing the loss effect in the
84
metal [86]. Here LD is the SPP dispersion length and β1 is the inverse of the SPP
wavevector.
For weak absorption where Γ ≪ 1, it is common to use perturbation method
to analytically approximate the result [101]. For the case of an SPP pulse on a
silver surface, however, the loss term can not be always regarded weak and the
numerical solution is the only way to demonstrate the evolution of the SPP pulse.
The numerical recipe of Crank Nicolson (CN) is used and for details one can refer
to Appendix D for the numerical algorithm and Appendix E for the source code
written in MatLab®.
To account for the metal loss more precisely, we use here a modified Drude
model for silver:
εm(ω) = ε∞ −
ω2p
ω2 + iγpω
+ i
σ
ε0ω
(V.56)
where γp is the damping coefficient and σ is the AC conductivity constant at optical
frequency. The physical meanings of all symbols and their values for silver are
summarized in table Tab.D.2 in Appendix D.
The SPP pulse duration is chosen to be 8fs so that the SPP linear dispersion
effect, represented by LD, and the nonlinear effect, represented by LNL, are
comparable to each other, both about 700µm. The repetition rate is tuned down
to 1Hz to avoid inadvertent heating caused by high peak intensities, in possible
experimental setups. The distance in the simulation is x = 2LD and when the
initial pulse is set to be exactly the soliton threshold for the lossless case, the pulse
85
Figure 5.7 The evolution of a dark SPP pulse with loss. The dark SPP pulse
travels on a metal surface with linear loss. Starting with a soliton profile, which
is a stable solution in lossless case, the pulse decreases smoothly with no internal
oscillations or break-ups.
is damped while propagating, but keeps a single pulse as shown in Fig. 5.7. With
other input amplitude distribution, a pulse break-up or internal oscillation occurs,
as shown in Fig. 5.8, due to the multiple breakup of high-order solitons and soliton
modulation instability [86] [102]. When the initial pulse amplitude is greater than
the soliton threshold, obvious pulse break-up occurs, similar to the lossless case in
Fig. 5.6a. After a short propagation distance, however, the amplitude of the pulse
86
decreases due to the loss, and the pulse break-up gradually turns into oscillations
as shown in Fig. 5.8a. When the initial pulse amplitude is less than the soliton
threshold, small oscillations propagate on both sides of the main pulse as shown in
Fig. 5.8b, similar to the lossless case in Fig. 5.6b.
(a) (b)
Figure 5.8 A schematic of pulse propagation with loss when the initial soliton
condition is not satisfied. (a) The optical pulse breaks up when the initial amplitude
is greater than the threshold. The pulse break-up is diminished during propagation
and gradually turns into oscillations because of the amplitude decrease due to the
loss (b) The optical pulse oscillates when the initial amplitude is less than the
threshold.
The damping of the dark SPP pulse may be compensated to keep a stable
soliton, by adding optical gains in the dielectric layer. Another possible way to
87
avoid the loss is by coherently coupling the dark SPP pulse to other plasmon modes,
with an effect similar to Electromagnetically induced transparency (EIT) [103].
Conclusions
In this chapter, we achieved a nonlinear Schro¨dinger (NLS) equation for a
temporal SPP pulse. The nonlinear dispersion relation of the temporal SPP pulse
was derived. This SPP pulse was supported by a nonlinear metal and a linear
dielectric. In the analytical studies, the nonlinear metal was dispersive and lossless,
while the linear dielectric was chosen to be vacuum. An analytical solution of the
NLS equation was achieved and was a temporal dark soliton for the nonlinear SPP
pulse. The soliton was infinitely stable during propagation for the ideal case. To
form this stable dark soliton, the relations among pulse peak amplitude, pulse
amplitude outline, pulse duration and background amplitude were achieved. For
a realistic metal with finite loss, numerical calculations were performed. With the
initial pulse satisfying the lossless soliton condition, the pulse maintained a single
recognizable peak when the amplitude was damped during propagation. With the
initial pulse not satisfying the lossless soliton condition, the pulse suffered internal
oscillation or pulse break-up for intensities too low or too high.
88
CHAPTER VI
CONCLUSIONS
The study of plasmonics has been drawing increased attention in recent
years. The physics of surface plasmon polaritons offers new and efficient ways
to guide light signals in nano-scale dimensions and at terahertz operation speeds
using metallic nanostructures. It also brings new technologies, such as chemical
sensing and biological sensing, to the other fields of science. To fully explore the
potential advantages of SPPs with respect to fast operation speed, broad bandwidth
and nanometer scale size, fundamental research into the nonlinear properties of
SPPs has become an active field, with its own challenges and excitements. This
dissertation serves to reveal a small tip of the iceberg.
We started this dissertation with a review of the history and recent
developments of plasmonics. In Chapter II, we presented a complete and detailed
derivation of major properties of linear surface plasmon polaritons (SPPs), such
as the dispersion relations. This derivation is the starting point of the following
chapters and it lays the foundation for our research. The SPPs were shown as a
guided wave between a dielectric and a metal surface. The optical properties of
the metals were analytically described by the Drude model when the frequencies of
the light waves were below the plasma frequency. A simple layout of planar SPPs
89
was used and the vacuum was used for the dielectric layer. The SPP dispersion
relations and the key dimensions of the SPPs were achieved.
In Chapter III, we introduced nonlinear optics in general then applied the
nonlinear wave equations to our specific nonlinear SPP system. The noble metals
in the SPP system possess third order nonlinearity, denoted with the third order
susceptibility χ(3). This χ(3) lead to the nonlinear polarization term PNL and then a
nonlinear wave equation for the metals was derived. In order to solve the nonlinear
wave equation, two specific cases were considered: sum frequency generation and
difference frequency generation. In the sum frequency generation case, the non-
depletion approximation was used to decouple the tripled frequency components
from the original frequency components. In the difference frequency generation
case, a nonlinear diffractive index was introduced to take into account the fact that
the generated frequency was indistinguishable from the original source frequency.
For both cases, a designated orientation of the metal crystal lattice was chosen to
reduce the nonlinear response of χ(3) to a scalar.
In Chapter IV, the steady-state solution of the nonlinear wave equation with
sum frequency generation was achieved. The solution showed third harmonic
generation in the nonlinear SPP process. The induced tripled frequency signal
from the SPPs is still bounded in the dielectric side of the interface. When the
optical properties of the metal are described by the classical Drude model, this
tripled frequency signal is able to propagate inside the metal with an obvious angle
90
of about 40◦ relative to the interface. The energy conversion efficiency from the
SPP to the tripled signal was analytically calculated. The efficiency is about 0.1%
using a standard continuous wave laser with visible light intensity I = 103W/cm2.
This low conversion efficiency complies with the non-depletion approximation made
when solving the nonlinear wave equations. The propagation angle of the tripled
signal was also calculated and found to have a dependency on the permittivity of
the metal and the incident wavelength of the SPP, and the signal is well separated
in space by its propagation direction from the original SPP waves. The angle ranges
from about 46◦ to about 56◦ when the incident wavelength is in the optical range.
For realistic metals, such as noble metals with interband transitions, absorption
above the plasma frequencies had to be considered. Tabulated data in literature
shows a transparency window around the plasma frequency for silver, with reduced
absorption, and consequently we used this transparency data to show possible
improvements for experimental setups, e.g. we considered a silver film with finite
thickness instead of a semi-infinite bulk. Within this transparency window, the
tripled signal has a transmittance of about 10% for a 50nm silver film, and the
propagation angle is about 40◦.
In Chapter V, we investigated the temporal properties of nonlinear SPPs
by considering the nonlinear Kerr effect in silver. First we derived the nonlinear
dispersion relation from the boundary conditions involving metal nonlinearities.
With the nonlinear dispersion relation, we applied slowly the varying amplitude
91
approximation and perturbed the wavevectors in the wave equations, and a
temporal nonlinear Schro¨dinger (NLS) equation for the nonlinear SPP propagation
was achieved. Given the values of the coefficients in the NLS equation, we found
that the stable solution is a SPP pulse, and that the shape of the pulse possesses
the form of a dark temporal soliton. The ideal dark soliton solution, when the
metal is lossless, was also demonstrated by a numerical calculation. When losses in
the metal accounted for, the numerical calculation showed that the evolution of the
pulse depends strongly on its shape. When the initial shape of the pulse satisfies
the same conditions of a soliton in the lossless case, the pulse propagated with only
broadening and damping. Otherwise, the pulse suffered background oscillations or
multi-peak breakups within a short distance after being launched.
The topics discussed in this dissertation are nonlinear plasmonics with the
consideration of the metallic nonlinearities. Nonlinear SPPs have many more
interesting properties to be explored. For future studies, we propose several
improvements and directions of efforts. For the THG of nonlinear SPPs, the energy
conversion efficiency can be increased by an incident SPP with higher intensity. The
power draining effect of the source SPP has to be considered when the efficiency
goes up where the non-depletion approximation is no longer valid, and the coupled
wave equations with both fundamental and tripled frequencies have to be solved
numerically. Without the non-depletion approximation, the propagation angle of
the tripled frequency signal can have dependency on the input optical intensity. In
92
addition, both the conversion efficiency and the propagation angle of the tripled
signal can be tunable if the dielectric layer is optically nonlinear or possesses
optical gain. Different structures are worth attention too. The proposed dielectric-
metal-dielectric (DMD) and metal-dielectric-metal (MDM) [104][105][106] layered
structures, in which each layer may be linear or nonlinear, may add exciting new
nonlinear effects in addition to what we have discussed in this dissertation. For
instance, the dispersion relation of the tripled signal in the dielectric side within
a single metal-dielectric interface may be altered by a nonlinear DMD or MDM
structure due to the coherence and interference of the optical waves in different
layers.
As for the SPP solitons, the main problem for achieving stable and long
distance SPP solitons is material loss. Loss compensation has been a major
direction of study in recent years, mostly in linear plasmonics [107]. A gain
dielectric medium may be capable of compensating for the metal loss. When the
dielectric layer has gain, the nonlinear boundary conditions are changed and thus
leads to a revised dispersion relation of the SPP. It is possible, with a gain medium
compensating the metal loss, that the nonlinear SPP solitons are stable and can
propagate long distance. It may even be possible that the dark soliton solution is
converted to a bright soliton solution, such as when the overall GVD is changed
from normal to anomalous, by the gain dielectric medium.
Nonlinear plasmonics provide us powerful tools to manipulate and guide
93
optical signals. For instance, devices with nanometer scale and terahertz operation
speed, as well as sub-wavelength digital signal transmission, may be utilized by
nonlinear SPP processes discussed in this dissertation. For the full potential of
the nonlinear plasmonics, it is necessary to carry out further systematic studies of
the nonlinear plasmonics, in both theory and experiment. The future studies of
the nonlinear plasmonics include, but not limited to, multi-layer structures such as
MDM and DMD, patterned guiding surfaces and supporting dielectric materials
with different optical properties such as optical nonlinearity and gain. Research on
nonlinear plasmonics not only improves understanding of the fundamental nature
of SPPs, it will also contribute novel methods to the fast growing technology world.
94
APPENDIX A
SOLUTION OF INHOMOGENEOUS SECOND ORDER ORDINARY
DIFFERENTIAL EQUATIONS
In the section, we solve an inhomogeneous second order ordinary differential
equation [108].
A second-order (inhomogeneous) ordinary linear differential equation of
function q(τ) encountered often has the form
d2q
dτ 2
+ 2ζ
dq
dτ
+ ω20q = Be
iωsτ (A.1)
where B is an arbitrary constant and ωs is the driving frequency.
The characteristic equation reads
r2 + 2ζr + ω20 = 0 (A.2)
The full solution of (A.1) comprises two parts: q = qt + qs where qt is the
general solution of the homogeneous equation, a.k.a the transient solution, while qs
is a particular solution, a.k.a. the steady-state solution.
For small damping |ζ | ≪ 1, the discriminant ∆ = (2ζ)2 − 4ω20 < 0, so the
characteristic equation has two solutions that are complex conjugate of one another:
r1 = −ζ + i
√
ω20 − ζ2 (A.3a)
r2 = −ζ − i
√
ω20 − ζ2 (A.3b)
95
Then the general (transient) solution to the homogeneous differential equation
is
qt(τ) = e
−ζτ
(
C1e
−iτ
√
ω2
0
−ζ2 + C2e
iτ
√
ω2
0
−ζ2
)
(A.4)
where C1 and C2 are determined by initial conditions. Note that qt → 0 when
t→∞ no matter how small the damping ζ is.
The particular (steady-state) solution qs can be found by the method of
undetermined coefficients. By making ansatz qs(τ) = A e
iωsτ where A is a complex
amplitude to be determined. By substituting qs back to (A.1) one finds A =
B/(−ω2s + 2iζωs + ω20). The particular solution thus is oscillatory
qs(τ) =
Beiωsτ
−ω2s + 2iζωs + ω20
(A.5)
In a real experiment, the damping coefficient is always finite, which gives
q(τ)
τ→∞
=⇒ qs(τ) (A.6)
Therefore, in the long time limit, only the steady-state solution is detected since
the general solution dies out soon.
96
APPENDIX B
HELMHOLTZ THEOREM AND HELMHOLTZ EQUATION
The Helmholtz theorem states [109]:
Any sufficiently smooth, rapidly decaying vector field can be resolved into
irrotational (curl-free) and solenoidal (divergence-free) component vector fields.
Given such a vector field V(r) satisfying
[∇ ·V]∞ = 0 (B.1a)
[∇×V]∞ = 0 (B.1b)
then there exist a scalar potential φ and a vector potential A such that
V(r) = −∇φ+∇×A
≡ V‖(r) +V⊥(r) (B.2)
where V‖ is the longitudinal (irrotational, curl-free) part that
∇×V‖ = 0 (B.3)
and V⊥ is the transverse (solenoidal, divergence-free) part that
∇ ·V⊥ = 0 (B.4)
Applying Fourier transformation on vector field V from r to k:
V˜(k) = F [V(r)] =
1(√
2pi
)3 ∫∫∫
r
V(r)eik·rd3r (B.5)
97
gives
k× V˜‖(k) = 0 (B.6a)
k · V˜⊥(k) = 0 (B.6b)
and consequently
∇×∇×V(r) =⇒
ik×
[
ik× V˜(k)
]
=ik×
[
ik× V˜⊥(k)
]
=ik
[
ik · V˜⊥(k)
]
− (ik)2V˜⊥(k)
=− (ik)2V˜⊥(k) (B.7)
98
APPENDIX C
THE ROTATION MATRIX
The nonlinear polarization is induced by the incident fields via the nonlinear
susceptibility χ(3) tensor as in (III.18). The first three terms on the right hand of
(III.18) are invariant under rotation thus independent of spatial coordinates and
the sum is labeled as PI where the superscript indicates “isotropic”. The last term
on the right hand of (III.18) is where the tensor nature of the χ(3) imbedded. It is
labeled as PA indicating “anisotropy” which is coordinate specific and defined as
PAi ≡ (M :: EaEbEc)i =
∑
jkl
M ijklE
a
jE
b
kE
c
l (C.1)
y
x
ϕ
ϕ
β
α
Eα = Ex cosϕ+ Ey sinϕ
Eβ = −Ex sinϕ+ Ey cosϕ
and
PAx = P
A
α cosϕ− PAβ sinϕ
PAy = P
A
α sinϕ+ P
A
β cosϕ
Figure C.1 The rotation matrix for (simplified 2-D) crystal - lab coordinates
conversion. The conversion relations for E and P are explicitly represented.
99
where ijkl represent any set of orthogonal coordinate axes. Now the nonlinear
polarization is separated to isotropic and anisotropic parts:
PNL = PI +PA (C.2)
Practically it is convenient to write both P and E in lab reference frame, so the
explicit form of tensor M in the lab coordinate frame is desired, and our central
goal of this section is to find the expressions for the tensor M .
Firstly we deal with the rotation part in the tensor M and add the effect of χ(3)
later. It is easier to break the total rotation matrix K between the cubic and the
lab frames to two steps: the rotation matrix T between the cubic and the crystal
frames, and the rotation matrix R between the crystal and the lab frames:
Ex′y′z′ = RExyz (C.3a)
Eαβγ = T
−1Ex′y′z′ (C.3b)
Px′y′z′ = TPαβγ (C.3c)
Pxyz = R
−1Px′y′z′ (C.3d)
Here E is the electric field vector, Exyz is the vector E expressed in lab frame, Ex′y′z′
is the vector E expressed in the crystal frame, and Eαβγ is the vector E expressed
in the cubic frame. The the rotation matrix R between lab and crystal frames with
rotation angle ϕ is defined as
R =
 cosϕ sinϕ 0− sinϕ cosϕ 0
0 0 1
 (C.4)
100
while the representation of rotation matrix T between the cubic and the crystal
frames depends on the orientation of cubic symmetry axes:
T(111) =

√
2
3
− 1√
6
− 1√
6
0
1√
2
− 1√
2
1√
3
1√
3
1√
3
 (C.5a)
T(110) =

− 1√
2
1√
2
0
0 0 1
1√
2
1√
2
0
 (C.5b)
T(001) =
1 0 00 1 0
0 0 1
 (C.5c)
The subscripts here indicate the orientation of the crystal, e.g. the direction of
the crystal z′ axis, while the actual crystal “cut” or “shearing plane” is in a plane
normal to this chosen axis, indicated as the shaded planes in Fig. C.4 Fig. C.3 and
Fig. C.2. The T is derived by simple geometric relations, for example the (111)
orientation gives:
xˆ′ =
√
2
3
xˆ− 1√
6
(yˆ + zˆ)
yˆ′ =
1√
2
(yˆ − zˆ)
zˆ′ =
1√
3
(xˆ+ yˆ + zˆ)
(C.6)
where we have chosen this definition so that the new x axis is projected on to the
original crystal axis in the plane of the crystal surface, shown in Fig. C.2.
101
xz
y
x
′
y′
z′
 = 1√
6
 2 −1 −10 √3 −√3√
2
√
2
√
2

xy
z

Figure C.2 The relation between (111) crystal “cut” (shaded) and the lab
coordinates. The direction of the crystal orientation, z′, is defined as the normal
to the shaded surface, while the other axes x′ and y′ are chosen with a degree of
freedom as long as the three axes comprise a right handed coordinate system. The
cubic coordinates are not shown here.
For the (110) surface we have
xˆ′ =
1√
2
(yˆ − xˆ)
yˆ′ = zˆ
zˆ′ =
1√
2
(xˆ+ yˆ)
(C.7)
and the matrix form is expressed in Fig. C.3. And finally for the (001) orientation
we simply choose the z′ axis to align normal to the surface which coincides with the
z axis, while x′ and y′ coincide with x and y in the lab frame. To directly connect
the cubic and lab coordinates, the compound matrix K ≡ R−1T is defined:
K(001) =
cosϕ − sinϕ 0sinϕ cosφ 0
0 0 1
 (C.8a)
102
z′
x′
y′
x
z
y
x
′
y′
z′
 =

− 1√
2
1√
2
0
0 0 1
1√
2
1√
2
0

xy
z

Figure C.3 The relation between (110) crystal “cut” (shaded) and the lab
coordinates. The direction of the crystal orientation, z′, is defined as the normal to
the shaded surface, shown in dotted line as well as x′ axis, while y′ points straight
up. The cubic coordinates are not show here.
K(111) =

√
2
3
cosϕ −cosϕ√
6
− sinϕ√
2
−cosϕ√
6
+
sinϕ√
2√
2
3
sinϕ −sinϕ√
6
+
cosϕ√
2
−sinϕ√
6
− cosϕ√
2
1√
3
1√
3
1√
3
 (C.8b)
K(110) =

−cosϕ√
2
cosϕ√
2
− sinϕ
−sinϕ√
2
sinϕ√
2
− cosϕ
1√
2
1√
2
0
 (C.8c)
The relation in the lab coordinate systems between the E and the P fields are
achieved by three steps:
1) From Exyz to Eαβγ by the rotation matrix:
Eαβγ = K
−1Exyz (C.9)
103
xz
y
x
′
y′
z′
 =
1 0 00 1 0
0 0 1

xy
z

Figure C.4 The relation between (100) crystal “cut” (shaded) and the lab
coordinates. The direction of the crystal orientation coincides with the axes of
the lab frame thus renders an identity matrix.
2) From Eαβγ to Pαβγ by χ
(3).
3) From Pαβγ to Pxyz by the rotation matrix:
Pxyz = KPαβγ (C.10)
and it is used to derived the explicit form of tensor M .
PAαβγ ≡
P
A
α
PAβ
PAγ
 =
E
a
αE
b
αE
c
α
EaβE
b
βE
c
β
EaβE
b
βE
c
β
 =
(K
′Eaxyz)α(K
′Ebxyz)α(K
′Ecxyz)α
(K ′Eaxyz)β(K
′Ebxyz)β(K
′Ecxyz)β
(K ′Eaxyz)γ(K
′Ebxyz)γ(K
′Ecxyz)γ
 (C.11)
and then in lab coordinates,P
A
x
PAy
PAz
 ≡ PAxyz = KPAαβγ = K
P
A
α
PAβ
PAγ

= K
(K
′
ααE
a
x +K
′
αβE
a
y +K
′
αγE
a
z )(K
′
ααE
b
x +K
′
αβE
b
y +K
′
αγE
b
z)(K
′
ααE
c
x +K
′
αβE
c
y +K
′
αγE
c
z)
(K ′βαE
a
x +K
′
ββE
a
y +K
′
βγE
a
z )(K
′
βαE
b
x +K
′
ββE
b
y +K
′
βγE
b
z)(K
′
βαE
c
x +K
′
ββE
c
y +K
′
βγE
c
z)
(K ′γαE
a
x +K
′
γβE
a
y +K
′
γγE
a
z )(K
′
γαE
b
x +K
′
γβE
b
y +K
′
γγE
b
z)(K
′
γαE
c
x +K
′
γβE
c
y +K
′
γγE
c
z)

104
= K
(KααE
a
x +KβαE
a
y +KγαE
a
z )(KααE
b
x +KβαE
b
y +KγαE
b
z)(KααE
c
x +KβαE
c
y +KγαE
c
z)
(KαβE
a
x +KββE
a
y +KγβE
a
z )(KαβE
b
x +KββE
b
y +KγβE
b
z)(KαβE
c
x +KββE
c
y +KγβE
c
z)
(KαγE
a
x +KβγE
a
y +KγγE
a
z )(KαγE
b
x +KβγE
b
y +KγγE
b
z)(KαγE
c
x +KβγE
c
y +KγγE
c
z)

(C.12)
where again PAxyz is the nonlinear polarization in lab frame and P
A
αβγ is the
nonlinear polarization in the cubic frame, with the superscript A representing the
anisotropic part of the polarization only. So we see, for instance we can determine
the elements of Mxxxx,M
x
yyy,M
x
zzz, · · · ,Mxxyz, · · · ,Mzxyz by
PAx =KααP
A
α +KαβP
A
β +KαγP
A
γ
=(K4αα +K
4
αβ +K
4
αγ)EEE|abcxxx
+ (KααK
3
βα +KαβK
3
ββ +KαγK
3
βγ)EEE|abcyyy
+ (KααK
3
γα +KαβK
3
γβ +KαγK
3
γγ)EEE|abczzz
+ (K3ααKβα +K
3
αβKββ +K
3
αγKβγ)EEE|abcxxy,xyx,yxx
+ (K3ααKγα +K
3
αβKγβ +K
3
αγKγγ)EEE|abcxxz,xzx,zxx
+ (K2ααK
2
βα +K
2
αβK
2
ββ +K
2
αγK
2
βγ)EEE|abcxyy,yxy,yyx
+ (K2ααK
2
γα +K
2
αβK
2
γβ +K
2
αγK
2
γγ)EEE|abcxzz,zxz,zzx
+ (KααK
2
βαKγα +KαβK
2
ββKγβ +KαγK
2
βγKγγ)EEE|abcyyz,yzy,zyy
+ (KααKβαK
2
γα +KαβKββK
2
γβ +KαγKβγK
2
γγ)EEE|abcyzz,zyz,zzy
+ (K2ααKβαKγα +K
2
αβKββKγβ +K
2
αγKβγKγγ)EEE|abcxyz,xzy,yxz,yzx,zxy,zyx
(C.13)
105
where shorthand notations are used for the summations of terms with permutation
subscripts, for example
EEE|abcxxz,xzx,zxx ≡ EaxEbxEcz + EaxEbzEcx + EazEbxEcx (C.14)
The other two components of PA are written in the same way as
PAy =KβαP
A
α +KββP
A
β +KβγP
A
γ
=(KβαK
3
αα +KββK
3
αβ +KβγK
3
αγ)EEE
abc
xxx
+ (K4βα +K
4
ββ +K
4
βγ)EEE|abcyyy
+ (KβαK
3
γα +KββK
3
γβ +KβγK
3
γγ)EEE|abczzz
+ (K2ααK
2
βα +K
2
αβK
2
ββ +K
2
αγK
2
βγ)EEE|abcxxy,xyx,yxx
+ (KβαK
2
ααKγα +KββK
2
αβKγβ +KβγK
2
αγKγγ)EEE|abcxxz,xzx,zxx
+ (KααK
3
βα +KαβK
3
ββ +KαγK
3
βγ)EEE|abcxyy,yxy,yyx
+ (KβαKααK
2
γα +KαβKββK
2
γβ +KβγKαγK
2
γγ)EEE|abcxzz,zxz,zzx
+ (K3βαKγα +K
3
ββKγβ +K
3
βγKγγ)EEE|abcyyz,yzy,zyy
+ (K2βαK
2
γα +K
2
ββK
2
γβ +K
2
βγK
2
γγ)EEE|abcyzz,zyz,zzy
+ (KααK
2
βαKγα +KαβK
2
ββKγβ +KαγK
2
βγKγγ)EEE|abcxyz,xzy,yxz,yzx,zxy,zyx
(C.15)
PAz =KγαP
A
α +KγβP
A
β +KγγP
A
γ
=(KγαK
3
αα +KγβK
3
αβ +KγγK
3
αγ)EEE|abcxxx
106
+ (KγαK
3
βα +KγβK
3
ββ +KγγK
3
βγ)EEE|abcyyy
+ (K4γα +K
4
γβ +K
4
γγ)EEE|abczzz
+ (KγαK
2
ααKβα +KγβK
2
αβKββ +KγγK
2
αγKβγ)EEE|abcxxy,xyx,yxx
+ (K2ααK
2
γα +K
2
αβK
2
γβ +K
2
αγK
2
γγ)EEE|abcxxz,xzx,zxx
+ (KγαKααK
2
βα +KγβKαβK
2
ββ +KγγKαγK
2
βγ)EEE|abcxyy,yxy,yyx
+ (KααK
3
γα +KαβK
3
γβ +KαγK
3
γγ)EEE|abcxzz,zxz,zzx
+ (K2βαK
2
γα +K
2
ββK
2
γβ +K
2
βγK
2
γγ)EEE|abcyyz,yzy,zyy
+ (KβαK
3
γα +KββK
3
γβ +KβγK
3
γγ)EEE|abcyzz,zyz,zzy
+ (KααKβαK
2
γα +KαβKββK
2
γβ +KαγKβγK
2
γγ)EEE|abcxyz,xzy,yxz,yzx,zxy,zyx
(C.16)
Now as a simple example for (100) symmetry axes, we can find the relation for
the representations of E and P fields the lab frame. Let us firstly evaluate PAx (the
shorthand notation for products of fields is used, such that Exyy ≡ ExEyEy):
PAx = P
A
α cosϕ− PAβ sinϕ
= EaαE
b
αE
c
α cosϕ−EaβEbβEcβ sinϕ
= (Eax cosϕ+ E
a
y sinϕ)(E
b
x cosϕ + E
b
y sinϕ)(E
c
x cosϕ+ E
c
y sinϕ) cosϕ
− (−Eax sinϕ+ Eay cosϕ)(−Ebx sinϕ+ Eby cosϕ)(−Ecx sinϕ+ Ecy cosϕ) sinϕ
= EaxE
b
xE
c
x cos
4 ϕ+ EaxE
b
xE
c
y cos
3 ϕ sinϕ+ EaxE
b
yE
c
x cos
3 ϕ sinϕ+ EaxE
b
yE
c
y cos
2 ϕ sin2 ϕ
+EayE
b
xE
c
x cos
3 ϕ sinϕ+ EayE
b
xE
c
y cos
2 ϕ sin2 ϕ+ EayE
b
yE
c
x cos
2 ϕ sin2 ϕ + EayE
b
yE
c
y cosϕ sin
3 ϕ
107
− [−EaxEbxEcx sin4 ϕ+ EaxEbxEcy cosϕ sin3 ϕ+ EaxEbyEcx cosϕ sin3 ϕ− EaxEbyEcy cos2 ϕ sin2 ϕ
+EayE
b
xE
c
x cosϕ sin
3 ϕ− EayEbxEcy cos2 ϕ sin2 ϕ− EayEbyEcx cos2 ϕ sin2 ϕ
+EayE
b
yE
c
y cos
3 ϕ sinϕ]
= Eabcxxx(cos
4 ϕ+ sin4 ϕ) + Eabcxxy(cos
3 ϕ sinϕ− sin3 ϕ cosϕ) + Eabcxyx(cos3 ϕ sinϕ− cosϕ sin3 ϕ)
+Eabcxyy(cos
2 ϕ sin2 ϕ+ cos2 ϕ sin2 ϕ)
+ Eabcyxx(cos
3 ϕ+ cosϕ sin3 ϕ) + Eabcyxy(cos
2 ϕ sin2 ϕ+ sin2 ϕ cos2 ϕ)
+Eabcyyx(cos
2 ϕ sin2 ϕ− cos2 ϕ sin2 ϕ)
+Eabcyyy(cosϕ sin
3 ϕ− cos3 ϕ sinϕ) (C.17)
then we get the matrix elements for M in (100) crystal face,
x, xxx(PAx ← Exxx) : cos4 ϕ+ sin4 ϕ =
1
4
(3 + cos 4ϕ) (C.18)
x, yyy(PAx ← Eyyy) : cosϕ sin3 ϕ− cos3 ϕ sinϕ = −
1
4
sin 4ϕ (C.19)
and they agree with the expression derived from the general discussion:
Mxxxx = K
4
αα +K
4
αβ +K
4
αγ = cos
4 ϕ+ sin4 ϕ =
1
4
(3 + cos 4ϕ) (C.20)
Mxyyy = KααK
3
βα +KαβK
3
ββ +KαγK
3
βγ = cosϕ sin
3 ϕ+ (− sinϕ) cos3 ϕ+ 0 = −1
4
sin 4ϕ
(C.21)
Similarly, as another example, for (111) crystal symmetry axes,
Mxxxz = K
3
ααKγα +K
3
αβKγβ +K
3
αγKγγ
=
(√
2
3
cosϕ
)3
1√
3
+
(
−cosϕ√
6
− sinϕ√
2
)3
1√
3
+
(
−cosϕ√
6
+
sinϕ√
2
)3
1√
3
108
=
1√
3
[
2
3
√
2
3
cos3 ϕ+ 2
(
−cosϕ√
6
)3
+ 6
(
−cosϕ√
6
)(
sinϕ√
2
)2]
=
1√
3
[
4√
6
cos3 ϕ− 3√
6
cosϕ
]
=
1√
18
cos 3ϕ (C.22)
And finally for (110) axes,
Mxxxx = K
4
αα +K
4
αβ +K
4
αγ
=
(
−cosϕ√
2
)4
+
(
cosϕ√
2
)4
+ (− sinϕ)4
=
1
2
cos4 ϕ+ sin4 ϕ
=
3
16
cos 4ϕ+
1
16
(9− 4 cos 2ϕ) (C.23)
Mzzzz = K
4
γα +K
4
γβ +K
4
γγ
=
(
1√
2
)4
+
(
1√
2
)4
+ 04
=
1
2
(C.24)
Myxxy = K
2
ααK
2
βα +K
2
αβK
2
ββ +K
2
αγK
2
βγ
=
(− cosϕ√
2
)2(− sinϕ√
2
)2
+
(
cosϕ√
2
)2(
sinϕ√
2
)2
+ (− sinϕ)2(cosϕ)2
=
3
16
(1− cos 4ϕ) (C.25)
Myxxx = KβαK
3
αα +KββK
3
αβ +KβγK
3
αγ
=
(
−sinϕ√
2
)(
−cosϕ√
2
)3
+
(
sinϕ√
2
)(
cosϕ√
2
)3
+ cosϕ(− sinϕ)3
=
sinϕ cos3 ϕ
2
− sin3 ϕ cosϕ
=
3
16
sin 4ϕ− 1
8
sin 2ϕ (C.26)
109
Finally the tensorM , representing the cubic rotation relative to the lab coordinates,
is determined for different forms for different cubic symmetry axe choices, i.e.
〈100〉, 〈111〉, 〈110〉
〈100〉 (C.27)
M ijkl xxx yyy zzz
x
3 + cos 4ϕ
4
−3 + sin 4ϕ
4
0
y
sin 4ϕ
4
3 + cos 4ϕ
4
0
z 0 0 1
(continued)
M ijkl xxy xxz xyy xzz yyz yzz xyz
x
sin 4ϕ
4
0
1− cos 4ϕ
4
0 0 0 0
y
1− cos 4ϕ
4
0 −sin 4ϕ
4
0 0 0 0
z 0 0 0 0 0 0 0
〈111〉 (C.28)
M ijkl xxx yyy zzz
x
1
2
0 0
y 0
1
2
0
z
cos 3ϕ√
18
−sin 3ϕ√
18
1
3
(continued)
M ijkl xxy xxz xyy xzz yyz yzz xyz
x 0
cos 3ϕ√
18
1
6
1
3
−cos 3ϕ√
18
0
sin 3ϕ√
18
y
1
6
sin 3ϕ√
18
0 0 −sin 3ϕ√
18
1
3
−cos 3ϕ√
18
z
sin 3ϕ√
18
1
3
−cos 3ϕ√
18
0
1
3
0 0
110
〈110〉 (C.29)
M ijkl xxx yyy zzz
x
9− cos 2ϕ+ 3 cos 4ϕ
16
−2 sin 2ϕ− 3 sin 4ϕ
16
0
y
−2 sin 2ϕ+ 3 sin 4ϕ
16
9 + 4 cos 2ϕ+ 3 cos 4ϕ
16
0
z 0 0
1
2
(continued)
M ijkl xxy xxz xyy
x
−2 sin 2ϕ+ 3 sin 4ϕ
16
0
3(1− cos 4ϕ)
16
y
3(1− cos 4ϕ)
16
0
−2 sin 2ϕ− 3 sin 4ϕ
16
z 0
1 + cos 2ϕ
4
0
(continued)
M ijkl xzz yyz yzz xyz
x
1 + cos 2ϕ
4
0
sin 2ϕ
4
0
y
sin 2ϕ
4
0
1− cos 2ϕ
4
0
z 0
1− cos 2ϕ
4
0
sin 2ϕ
4
For the special angle value ϕ = 0, M takes simple forms
〈100〉
M ijkl xxx yyy zzz xxy xxz xyy xzz yyz yzz xyz
x 1 0 0 0 0 0 0 0 0 0
y 0 1 0 0 0 0 0 0 0 0
z 0 0 1 0 0 0 0 0 0 0
(C.30)
111
〈111〉
M ijkl xxx yyy zzz xxy xxz xyy xzz yyz yzz xyz
x
1
2
0 0 0
1√
18
1
6
1
3
− 1√
18
0 0
y 0
1
2
0
1
6
0 0 0 0
1
3
− 1√
18
z
1√
18
0
1
3
0
1
3
− 1√
18
0
1
3
0 0
(C.31)
〈110〉
M ijkl xxx yyy zzz xxy xxz xyy xzz yyz yzz xyz
x
1
2
0 0 0 0 0
1
2
0 0 0
y 0 1 0 0 0 0 0 0 0 0
z 0 0
1
2
0
1
2
0 0 0 0 0
(C.32)
where the simplest case for (100) surface placed crystal with no rotation, the tensor
M takes the form:
M ijkl =

1 when jkl = xxx, yyy or zzz
0 otherwise
(C.33)
where i = x, y or z. So the anisotropic term in (III.18) is simplified to
χdM :: E
aEbEc = χd(xˆE
a
xE
b
xE
c
x + yˆE
a
yE
b
yE
c
y + zˆE
a
zE
b
zE
c
z) (C.34)
and we are going to use this simplest configuration in (C.33) or (C.34) to illustrate
the nonlinear SPP properties, without the cumbersome geometry distraction, in all
the related chapters.
112
APPENDIX D
NUMERICAL ALGORITHM FOR NLS EQUATIONS
The NLS equation with a loss term (V.55) is prototyped in the class
of convection-diffusion equations for the Finite difference method for partial
differential equations [110]. For both accuracy and efficiency in computation,
the implicit Crank-Nicolson scheme is adapted for the NLS equation [111]. For
notational convenience and comparison to heat equation, the NLS equation (V.55)
is rewritten with space-time variables switched, with first derivative on the left side
and other terms on the right side:
i
∂u
∂t
= −1
2
∂2u
∂x2
= |u|2u+ iΓu (D.1)
To discretize the equation, we use forward-time and central-space meshes and label
the time-mesh interval ∆t and the space-mesh interval ∆x. For example, the time
derivative term on the left side of (D.1) is
∂u
∂t
≃ u
n+1
k − unk
∆t
(D.2)
where the superscripts n + 1 and n label the time step with n(= 1, · · · , N) as the
current time step, and the subscript k(= 1, · · · , K) labels the central space point.
The last term on the right hand side of (D.1) is represented by its central space
113
point, but average of forward time steps:
u ≃ u
n+1
k + u
n
k
2
(D.3)
and similarly the second term:
|u|3 ≃ |u
n+1
k |3 + |unk |3
2
(D.4)
As for the second spatial derivative term
∂2u
∂x2
, we firstly write the standard central-
space derivative at current time step n:
∂2u
∂x2
≃ u
n
k+1 − unk + unk−1
∆x2
(D.5)
then replace every term un by its forward time average (un + un+1)/2:
∂2u
∂x2
≃ u
n
k+1 − unk + unk−1
2∆x2
+
un+1k+1 − un+1k + un+1k−1
2∆x2
(D.6)
Define the diffusion number, which is dimensionless and positive
λ =
∆t
∆x2
(D.7)
and collect all discrete expression for (D.1), we have now
2i(un+1k − unk) =−
λ
2
(un+1k+1 − un+1k + un+1k−1) + ∆t
(|un+1k |3 + iΓun+1k )
− λ
2
(unk+1 − unk + unk−1) + ∆t
(|unk |3 + iΓunk) (D.8)
Group the linear forward time steps n+1 to one side, while the current time steps
n and higher orders of n+ 1 to the other side:
λ
2
un+1k−1 +
[
2i−
(
λ
2
+ iΓ∆t
)]
un+1k +
λ
2
un+1k−1
=− λ
2
unk−1 +
[
2i+
(
λ
2
+ iΓ∆t
)]
unk −
λ
2
unk−1 +∆t
(|un+1k |3 + |unk |3)
(D.9)
114
To update all spatial points k = 1, · · · , K for time step n + 1, the previous time
step n values are used, as illustrated in Fig. D.1. Define 2Λ ≡ λ
2
+ iΓ∆t and use
x
t
n+ 1
Unknown
n
Known
n− 1 Known
k − 1 k k + 1
Figure D.1 Crank-Nicolson scheme illustrated with every point representing a
value of unk at time step n and position k. The circles are grid points that are
known at current time step and the squares are unknown points. The forward time
step point, indicated as a solid square (red), is determined by known points with
previous time step, indicated as shaded circles (blue), and unknown grid points,
indicated as hollow squares (red), at the same time step using an implicit method.
115
the matrix form:
2(i− Λ) λ
2
0 · · · 0 0 −λ
2
λ
2
2(i− Λ) λ
2
· · · 0 0 0
0
. . .
. . .
. . . 0 0
0 · · · λ
2
2(i− Λ) λ
2
· · · 0
0 0
. . .
. . .
. . . 0
0 0 · · · 0 λ
2
2(i− Λ) λ
2
−λ
2
0 · · · 0 0 λ
2
2(i− Λ)


u1
...
uk−1
uk
uk+1
...
uK

n+1
=

2(i+ Λ) −λ
2
0 · · · 0 0 λ
2
−λ
2
2(i+ Λ) −λ
2
· · · 0 0 0
0
. . .
. . .
. . . 0 0
0 · · · −λ
2
2(i+ Λ) −λ
2
· · · 0
0 0
. . .
. . .
. . . 0
0 0 · · · 0 −λ
2
2(i+ Λ) −λ
2
λ
2
0 · · · 0 0 −λ
2
2(i+ Λ)


u1
...
uk−1
uk
uk+1
...
uK

n
+∆t


|u1|3
...
|uk−1|3
|uk|3
|uk+1|3
...
|uK |3

n
+

|u1|3
...
|uk−1|3
|uk|3
|uk+1|3
...
|uK|3

n+1
(D.10)
where the common time step stamps n and n+1 are labeled out side of the column
vectors in the matrix expression. Accordingly we can define matrix A, column
116
vector U and matrix B to simplify the notation such that (D.10) is written as
AUn+1 = BUn +∆t
(|Un+1|3 + |Un|3) (D.11)
where the superscripts n and n + 1 label time steps while 3 means to the third
power. Now the Un+1 update formula is
Un+1 = A−1BUn + A−1∆t|Un|3 + A−1∆t|Un+1|3 (D.12)
The right side of the update scheme, however, contains a forward time term |Un+1|3.
Therefore, an implicit iteration method is used to determine the new time step value
for Un+1. Namely, to achieve the next time step Un+1, we start with a current time
step value Un for all U terms on the right side and get a tentative value Un+1(i = 0)
on the left side for the i = 0 iteration. Then the Un+1(i = 0) value is used back
on the right side to get an improved result for the left side Un+1(i = 1). Keep the
iteration to improve the forward time Un+1 term until i = ∞ or the accuracy is
within a small range of tolerance.
The expected solution is a dark soliton and it is a hole imbedded in a non-zero
background. Thus the boundary condition for the numerical method, instead
of commonly used periodic boundary condition [110] or transparent boundary
condition [112], an inverse-periodic boundary condition is used, which has a pi
phase shift between the left end and right end of the grid. This inverse-periodic
boundary condition is represented in the matrix A and B by the off-diagonal corner
117
Table D.1 Numerical solution with lossless metals.
c = 3× 108(m/s) Speed of light in vacuum
ε0 = 8.8542× 10−12(F/m) Vacuum permittivity
εd = 1 Dielectric constant of vacuum
ωp = 1.3814× 1016(rad/s) Plasma frequency of metal
ε∞ = 5.1 Background dielectric of metal
λ0 = 1500(nm) Laser wavelength
k0 =
2pi
λ0
= 4.18879× 106(rad/m) Laser wavevector in vacuum
ω0 =
2pic
λ0
= 1.25664× 1015(rad/s) Angular frequency of laser
β0 = β(ω0) = 4.21013× 106(rad/m) SPP wavevector(∼ kspp)
β1 = β
′(ω0) = 3.38631× 10−9(s/m) Inverse of SPP group velocity (1/vg)
β2 = β
′′(ω0) = 9.27626× 10−26(s2/m) SPP GVD
χ = χ(3) = 2.8× 10−19(m2/V 2) Nonlinear susceptibility of metal
ν =
k0χ
(3)
√
ε3d√
16|εm0|(|εm0| − εd)3
Nonlinear SPP wavevector variant
= 2.95× 10−17(m/V 2)
Area= 0.01(mm2) Laser spot size
τ0 = 10(ps) Laser pulse duration
L0 = τ0/β1 = 2.95× 10−3(m) Pulse length, 2.95mm
LD = τ
2
0 /β2 = 1078(m) SPP Dispersion length
frep = 80(MHz) Pulsed laser repetition value
A0 =
√
β2√
ντ0 cos φ
= 5.6× 106(V/m) Threshold amplitude for fundamental soliton
Pave = 0.338(W ) Average power (pulsed laser)
elements, which possess an additional negative sign comparing to their ”neighbor”
elements.
The parameters used in the numerical calculation and the results are tabulated
in Tab.D.1 and Tab.D.2
118
Table D.2 Numerical solution with lossy metals.
ωp = 1.3814× 1016(rad/s) Plasma frequency of silver
γp = 6.857× 1013(rad/s) Damping coefficient
ε∞ = 5.1 Background dielectric of metal
λ0 = 1500(nm) Laser wavelength
τ0 = 8(fs) Laser pulse duration
L0 = τ0/β1 = 2.36(µm) Pulse length
Lspp = 641(µm) SPP propagation length
LD = 741(µm) SPP Dispersion length
LNL =
1
Reν|A|2 = 679(µm) Nonlinear effect length
frep = 1(Hz) Pulsed laser repetition value
A0 =
√
β2√
ντ0 cos φ
= 7420(V/m) Threshold amplitude for fundamental soliton
Pave = 5.77(µW ) Pulse power
Γ = 2.31 Loss coefficient
%MatLab source code ;
%Source Code f o r Crank Nico l son Method − Crank Nico l son .m
c l e a r a l l ;
c l c ;
%c l f ;
%%
%Parameters
mesh_size = 1000;
total_step = 2000;
dx = 0 . 0 2 ;
dt = dx * dx / 5 . 0 ;
amplitude = 4 . 0 ;
%v e l o c i t y = 0 . 1 0 ;
grid_mem = 3 ;
Gamma = −0.0;%damping
%%
lambda=dt/dx/dx ;
119
diag_ele = lambda − 1i*dt*Gamma ;
off_ele = −0.5* lambda ;
% diagna l −1sub +1super
i=[1: mesh_size 2 : mesh_size 1 : mesh_size −1] ;
j=[1: mesh_size 1 : mesh_size−1 2 : mesh_size ] ;%( i , j )−>s
s=[diag_ele *ones (1 , mesh_size ) off_ele *ones (1 , mesh_size−1) . . .
off_ele *ones (1 , mesh_size−1) ] ;
imagi = 1i* speye ( mesh_size ) ;
major = spar se (i , j , s ) ;
leftA = imagi + major ;
leftA (1 , mesh_size ) = −off_ele ;
leftA ( mesh_size , 1 ) = −off_ele ;
rightB = imagi − major ;
rightB (1 , mesh_size ) = off_ele ;
rightB ( mesh_size , 1 ) = off_ele ;
c l e a r diag_ele off_ele i j s ;
invA=inv ( leftA ) ;
%%
%Graph Options
n_lines = 10;
line_skip = max(1 , uint32 ( total_step / n_lines ) ) ;
n_points = 100 ;
point_skip = max(1 , mesh_size / n_points ) ;
%%
%I n i t i a l Condit ion
x = (−mesh_size /2+1: mesh_size /2) *dx ;
init_line = amplitude * tanh ( x ) ;
%i n i t l i n e = −2 i * amplitude *exp(−2 i * v e l o c i t y *x) .* sech (2* amplitude *x ) ;
%%
%Grid prepare
grids = [ init_line ; init_line ; init_line ] ;
%
120
f o r time_step = grid_mem : 2* grid_mem
t3 = mod ( time_step − 1 , grid_mem ) + 1 ;%3
t2 = mod ( time_step − 2 , grid_mem ) + 1 ;%2
t1 = mod ( time_step − 3 , grid_mem ) + 1 ;%1
%FTCS
u2C = double ( grids (t2 , : ) ) ;%Center
u2L = double ( [ grids ( t2 , 1 ) grids ( t2 , 1 : end−1) ] ) ;%Lef t abrupt boundary cond i t i on
u2R = double ( [ grids ( t2 , 2 : end ) grids (t2 , end ) ] ) ;%Right
grids (t3 , : ) = double (1 i*lambda *(2* u2C−(u2L+u2R ) ) + . . .
(1+2i*dt*u2C .* conj ( u2C ) ) .* u2C ) ;
end
c l e a r t3 t2 t1 u2C u2L u2R
water = abs ( grids ( grid_mem , 1 : point_skip : mesh_size ) ) ;
%%
%evo l u t i on
f i g u r e (1) ;
c l f ;
p l o t (x , abs ( grids ( grid_mem , : ) ) )
hold on ;
f o r time_step = grid_mem : total_step
t3 = mod ( time_step − 1 , grid_mem ) + 1 ;%3 today
t2 = mod ( time_step − 2 , grid_mem ) + 1 ;%2 yes terday
t1 = mod ( time_step − 3 , grid_mem ) + 1 ;%1 the day be f o r e yes te rday
u2 = grids ( t2 , : ) ;
%Du Fort
% alpha = double (1 i /2/ dt + 1/dx/dx ) ;
% beta = double (1 i /2/ dt − 1/dx/dx ) ;
% u2L = double ( [ g r i d s ( t2 , 1 ) g r i d s ( t2 , 1 : end−1) ] ) ;%Lef t −
% u2R = double ( [ g r i d s ( t2 , 2 : end ) g r i d s ( t2 , end ) ] ) ;%Right
% u2C = gr i d s ( t2 , : ) ;
% u1C = double ( g r i d s ( t1 , : ) ) ;%Center
% gr i d s ( t3 , : ) = double ( ( beta *u1C + (u2L+u2R) /dx/dx − . . .
121
%2*u2C .* conj (u2C) .*u2C ) / alpha ) +1 i *Gamma*u2C/alpha ;
%
%Crank Nicol son
rightF1 = u2*rightB . ' ;
rightF2 = 0*dt .* u2 .* conj ( u2 ) .* u2 ;
rightFixed = rightF1+rightF2 ;
% vnewU=[u2 ; u2 ] ;% dimension r e s e r v a t i o n
% t o l e r = 1e−12;% to l e r ance
% i t e r o l d = 1 ;
% i t e r new = 2 ;
% e r r t = 1 ;
% whi le ( e r r t > t o l e r )
% oldU = vnewU( i t e r o l d , : ) ;
% rh s i d e = r ightF ixed + dt *oldU .* conj ( oldU ) .* oldU ;
% vnewU( i te r new , : ) = invA* rhs ide ' ;
% e r r t = max( abs ( (vnewU ( 1 , : )−vnewU ( 2 , : ) ) ) ) ;
% i t e r n ew = 3 − i t e r n ew ;
% i t e r o l d = 3 − i t e r o l d ;
% end
% gr i d s ( t3 , : ) = vnewU( i te r new , : ) ;
% f o r i t e r = 1 :3
rhside = rightFixed + 0*dt*u2 .* conj ( u2 ) .* u2 ;
u2 = rhside *invA . ' ;
% end
grids (t3 , : ) = u2 ;
i f ( mod ( time_step , line_skip ) == 0)
d i sp ( time_step ) ;
f i g u r e (1) ;
p l o t (x , abs ( grids (t3 , : ) ) ) ;
water = [ water ; abs ( grids ( t3 , 1 : point_skip : mesh_size ) ) ] ;
122
end
end
f i g u r e (2) ;
w a t e r f a l l ( water ) ;
hold off ;
%end
123
APPENDIX E
C/C++ SOURCE CODE FOR FDTD METHOD
FDTDmain.cpp
#include <cstdlib>
#include <iostream>
// This program is a simulation f o r SPP propagation on a metal s u r f a c e .
// The simulation f o r linear metal is stable
// The simulation f o r nonlinear metal is not stable f o r large nonlinear
// susceptibilities , because metal has negative dielectric constant
//which does not compatible with a Kerr nonlinear index .
// Reference Book : Allen Taflonve , Susan C . Hagness , Computational electrodynamics ,
//The f i n i t e −difference time−domain method , 3rd edition , page 203
//This program is f o r : Complex vs Real calculation optional .
//This program is ALL r e a l valued . Complex numbers are computed by separating
// r e a l and imaginary parts
//Unit : D=H*k/omega
//unit : SI ( modified )
#include ”math . h”
#include ”stdio . h”
#include ”stdlib . h”
#include ”algorithm ”
124
#include ”time . h”
#include ”FDTDarray . h”
#include ”FDTDcubic . h”
#include ”FDTDtime . h”
#define SQR2 1.414216562373// sq r t (2)
#define PI 3.1415926535898
using namespace std ;
// parameters definition , caution : only constants are defined as extern variables
/*
Naming convention notes :
1 . upper case f o r grids values
2 . ” _0” f o r constants
3 . ” _func ” f o r functions
4 . ” _r” f o r relative , dimensionless values
5 . prefix ”d” means difference
6 . prefix ”C” f o r factors
7 . postscript ”c” f o r complex values
*/ //
/*−−−−−−−−−−−−−−−−−−−−Universal Constants−−−−−−−−−−−−−−−−−−−−−*/
double const epsilon_0 =8.854187817 e−3;
//unit : nF/m , dielectric constant in vacuum . 1E−12−−>F/m
double const mu_0=0.4e3*PI ;
//unit : nH/m , permeability constant in vacuum . 1E−6−−>H/m
double const c_0=1/ sq r t ( epsilon_0 *mu_0 ) ;
//==0.29 , unit : nm/attosec , speed of light in vacuum . 1E−9/1E−18=1E9−−>m/s
double lambda_0=800;
//unit : nm , center optical wavelength . 1E−9−−>m
double const omega_0=2*PI*c_0/ lambda_0 ;
//unit : rad/attosec , central angular frequency , ( 1 E9*1E9=1E18−−>Hz ) ;
double const k_0=omega_0 /c_0 ;
//unit : /nm , wavevector . 1E9−−>/m
// typical optical (500 nm ) ˜ 2PI ( 0 . 0006 ) GGHz=2PI ( 6 . ( ( E−4) ) ) GGHz −>2PI ( 6 . E14 ) Hz
125
double const epsilon_inf =5.1;
// f o r silver , use relative values !
//−−−−−−−−−−For Debye Model ( Optional )−−−−−−−−−−−−−
double const epsilon_sp =1;
//−16643.2* epsilon_0 ; // pF/m . static response , −16643.2 f o r silver
// double const epsilon_inf =1;//5.1* epsilon_0 ; // infinite response , 10 .9866 f o r silver
double const depsilon_p =0;
// epsilon_sp−epsilon_inf ; // difference
double const tau_p=9076.707;
// relaxation time : attosec ; 9 . 076707 e−15 s f o r silver
//−−−−−−−−For Lorentz Model ( Optional )−−−−−−−−−−−
double const delta_p=0;
//1E−4*omega_p ; // damping coefficient , TBA
//−−−−−−−−−−−For Drude Model ( in use )−−−−−−−−−−−−−−
double const omega_p=1.3814E−2;
// frequency of pole pair , GGHz
double const gamma_p=3.18786E−6;
// relaxation time , unit GGHz
//−−−−−−−−Kerr Nonlinearity ( optional )−−−−−−−−
//ref . Nonlinear Optics by Rober W . Boyd , 3rd Edition , Page 222 ,212
// double const N_atom=4e22 ; / / ( cm )−3 atom number density
// double const e_charge=1.6e−19;// Coulomb , unit charge
// double const m_atom=1;
double const alpha=1;// alpha=1 f o r All Kerr and None Raman
double const chi3M=2.8e−19;//2.8e2 ; / / 2 . 8 ;
//unit ( m2/V2 )−−−−−Ag : 2 . 8 e−19 ( m2/V2 ) , Au 7 . 6 e−19
//−−−−−−−−Raman Effect ( Optional )−−−−−−−−
//ref . Nonlinear Optics by Rober W . Boyd , 3rd Edition , Page 222
double const a0_bohr =0;//0 . 05 ; // unit : nm , Bohr radius ,
double const v_electron =0;// c_0 /137;// nm/attosec , typical electronic velocity
double const tau1=0;
//2*PI*a0_bohr / v_electron ; //143 attosec , or 1E−16 second ;
126
//a single Lorentzian l i n e centered on the optical phonon frequency 1/t1 , TBA
double const tau2=0;// tau1 /1E4 ;
//and having a bandwidth of 1/t2 , the reciprocal phonon lifetime , TBA
double const omega_raman =0;// sq r t ( tau1*tau1+tau2*tau2 ) /tau1/tau2 ;
double const delta_raman =0;//1/tau2 ;
//−−−−−−−SPP parameters−−−−−−−−−−−−
double const drudedenom=(pow ( omega_0 , 4 )−pow ( omega_0 *gamma_p , 2 ) ) ;
double const epsilon_R=epsilon_inf−pow ( omega_p *omega_0 , 2 ) /drudedenom ;
// r e a l part of dielectric constant
double const epsilon_I=omega_0 *omega_p *omega_p *gamma_p /drudedenom ;
// imaginary part of dielectric constant
double const epsilon0R2=(1+epsilon_R ) *(1+ epsilon_R ) ; // squared
double const kx_R=k_0 *sqrtReal ( epsilon_R , epsilon_I ,1+ epsilon_R , epsilon_I ) ;
// wavevector , unit : 1/nm
double const kx_I=k_0 *sqrtImag ( epsilon_R , epsilon_I ,1+ epsilon_R , epsilon_I ) ;
double const ky0_R=k_0*sqrtReal (−1 ,0 ,1+ epsilon_R , epsilon_I ) ;
// decay constant , kappa , in air
double const ky0_I=k_0*sqrtImag (−1 ,0 ,1+ epsilon_R , epsilon_I ) ;
// oscillation in air
double const ky2_R=k_0*sqrtReal ( epsilon_I *epsilon_I−
epsilon_R *epsilon_R ,−2* epsilon_R *epsilon_I ,1+ epsilon_R , epsilon_I ) ;
// decay in metal
double const ky2_I=k_0*sqrtImag ( epsilon_I *epsilon_I−
epsilon_R *epsilon_R ,−2* epsilon_R *epsilon_I ,1+ epsilon_R , epsilon_I ) ;
// oscilation in metal
double const lambda=2*PI/kx_R ; // SPP wavelength along x ( propagation direction )
// double mu_r=1;// dimensionless , relative permeability constant to mu_0 ;
// double impedance_0 =376.730;// units : Ohm , s q r t ( epsilon_0 /mu_0 )
/*−−−−−−−−−−−−−−−−−−−−−−−−−FDTD parameters−−−−−−−−−−−−−−−−−−−−−−−−−−−−*/
int const Ntotal=100000;
127
//100000;// time steps number>3 41 ,141 ,241
int filestep=Ntotal /10 ;
//file written every ”filestep ” steps . Max=Ntotal , Min=1
int const halfgridsI =500;
//just f o r convienence , half incident area
int const gridsI=halfgridsI *2 ;
// simulation fields area
int const halfgridsJ =150;
// simulation fields area
int const gridsJ=halfgridsJ *2 ;
// simulation fields area−−−−−−−−−−−>x ( propagation direction )
int const CPML_n=30;
// Perfect Match Lay GRIDS width (30 GRIDS on each side )
int const TFSF_n=2;
//Total−Field/Scattered−Field technique layer width
int const SII=gridsI+CPML_n*2+TFSF_n *2 ;
// space GRIDS number , x ax i s
int const SJJ=gridsJ+CPML_n*2+TFSF_n *2 ;
// space GRIDS number , y ax i s
//int SKK=SII ; // space GRIDS nubmer , z ax i s ( only f o r 3D , not in use )
//z (H−field direction f o r TM mode , dependency not in use f o r 2D case )
int const centerI=SII /2 ;
int const centerJ=SJJ /2 ;
int const metalJ=SJJ /2;// metal−air boundary f o r y
int const SFt=SJJ−CPML_n ; // Top boundary f o r Scattered fields
int const SFb=CPML_n ; // Bottom , Scattered Fields
int const SFl=CPML_n ; // Left , Scattered fields
int const SFr=SII−CPML_n ; // Right boundary f o r Scattered fields
int const TFt=SFt ; // centerJ+halfgridsJ ; // top , Total fields
int const TFb=SFb ; // centerJ−halfgridsJ ; // bottom , T fields
int const TFl=SFl+TFSF_n ; // centerI−halfgridsI ; // left
int const TFr=SFr ; // centerI+halfgridsI ; // right
128
int const Fileoffset=0*TFSF_n+0*CPML_n ; // writing index
int const nonlinearM =100;// iteration loop number f o r nonlinear case 3 ;
double const ss_0=1/SQR2 ;
// stablity constant f o r 2D case , dt half index shall be considered
// double const omega_0=2*PI *1000E−6;
// double const lambda_0=2*PI*c_0/omega_0 ;
short Nlambday=50;// grids sampling density , at least 20 ,
short Nlambdax=200;// related to dt calculation , at least 50
double dymax=lambda_0 /Nlambday ;
double dxmax=lambda_0 /Nlambdax ;
double maxVSmin=lambda_0 *ky2_R ;
// ratio of dmax : dmin , wavelength / penetration in metal , ˜10
int varygrade_n =40;// number of vary nonuniform grids f o r metal
int mingrade_n =50;// number of minimun grids
double dxmin=dxmax ; //* pow ( gradeRatex , gradeMinLeft−gradeLeft ) ;
double dymin=dymax /maxVSmin ; //* pow ( gradeRatey , gradeMinDown−gradeDown ) ;
double gradeRatex =1;
double gradeRatey=pow ( maxVSmin , 1 . 0 / varygrade_n ) ;
//make sure adjacent cells width ratio between 0.5˜2
double const dt=1.0/(c_0 * s q r t (1/ dxmin/dxmin+1/dymin/dymin ) ) ;
//unit : atto second , use the smaller cell value , maximum dt value
double const phi=0;
//unit : radian , angle of incidence wave vector to x−axi s ,
//−90˜90 degree f o r this program
//******* FDTD extern parameters end here *******************/
// f o r units , c : 0 . 3 nm/attosecond , dt : attosec , dx : nm , f : 1 MTHz
double const a_debye=0;//(2* tau_p−dt ) /(2* tau_p+dt ) ;
double const b_debye=0;// epsilon_0 *depsilon_p *dt /(2* tau_p+dt ) ;
129
double const a_lorentz=0;//(2− omega_p *omega_p *dt*dt ) /( delta_p *dt+1) ;
double const b_lorentz =0;//( delta_p *dt−1) /( delta_p *dt+1) ;
double const c_lorentz=0;
//( epsilon_0 *depsilon_p *omega_p *omega_p *dt*dt ) /( delta_p *dt+1) ;
double const drudeFDTDdenom=2+dt*gamma_p ; // denominator in P [ n+1] calculation
double const a_drude=4/drudeFDTDdenom ;
double const b_drude=(dt*gamma_p −2)/drudeFDTDdenom ;
double const c_drude=epsilon_0 *2* omega_p *omega_p *dt*dt/drudeFDTDdenom ;
// epsilon_0 included !
double const a_raman=0;//(2− omega_raman *omega_raman *dt*dt ) /( delta_raman *dt+1) ;
double const b_raman=0;//( delta_raman *dt−1) /( delta_raman+1) ;
double const c_raman=0;
//((1− alpha ) *chi3* omega_raman * omega_raman *dt*dt ) /( delta_raman *dt+1) ;
int SFflag ( int ii , int jj ) ;
// Boolean f l a g . True 1 : ii , jj in scattered fields area . False 0 : not
int TFflag ( int ii , int jj ) ;
// Boolean f l a g . True 1 : ii , jj in total fields area . False 0 : not
int CPMLflag ( int ii , int jj ) ;
// Boolean f l a g . True 1 : ii , jj in CPML area . False 0 : not
int Metalflag ( int ii , int jj ) ;
int Airflag ( int ii , int jj ) ;
/** MAIN ********************************************************/
int main ( ) {// r eturn value=1: e r r o r . r e turn value=0: normal exit .
int nn=0;
int ii=0;
int jj=0;//ii , jj f o r space , nn f o r time
int pp=0;// f o r pole loops
int np0=0,np1=1,np2=2;// parity of nn , as time index f l a g
double nT=2*PI/omega_0 /dt ; // period of plane wave , in terms of time steps .
int FileFlag=0;
130
double ydis=0;// abs distance to orgin point x=0,y=0
// double xdis=0;
double ** Hzin=array2D (2 , SJJ+1) ; / / 1 . Incident wave source , ADE
double ** Eyin=array2D (2 , SJJ+1) ; / / 2 . Incident wave source , ADE
double *** Hz=array3D (3 , SII+1,SJJ+1) ;
//3 . space cell [ 2 ] [ x−g r i d index ] [ y−g r i d index ] f o r Hz which takes the edges
double *** Ey=array3D (3 , SII+1,SJJ+1) ; / / 4 . space cell [ 2 ] [ x ] [ y ] f o r Ey ,
double *** Ex=array3D (3 , SII+1,SJJ+1) ; / / 5 .
double *** Dy=array3D (3 , SII+1,SJJ+1) ; / / 6 . 2 Displacement y
double *** Dx=array3D (3 , SII+1,SJJ+1) ; / / 7 . 2 Displacement x
double *** Hzx=array3D (3 , SII+1,SJJ+1) ; / / 8 . Hz component
double *** Hzy=array3D (3 , SII+1,SJJ+1) ; / / 9 . Hz component
double *** Px_Debye=array3D (3 , SII+1,SJJ+1) ;//10 , polarization terms
double *** Px_Lorentz=array3D (3 , SII+1,SJJ+1) ;//11
double *** Px_Drude=array3D (3 , SII+1,SJJ+1) ;//12
double *** Py_Drude=array3D (3 , SII+1,SJJ+1) ;//13
double *** Px_Kerr=array3D (3 , SII+1,SJJ+1) ;//14
double *** Px_Raman=array3D (3 , SII+1,SJJ+1) ;//15
double *** Py_Debye=array3D (3 , SII+1,SJJ+1) ;//16
double *** Py_Lorentz=array3D (3 , SII+1,SJJ+1) ;//17
double *** Py_Kerr=array3D (3 , SII+1,SJJ+1) ;//18
double *** Py_Raman=array3D (3 , SII+1,SJJ+1) ;//19
double *** Sp=array3D (3 , SII+1,SJJ+1) ;//20
double ** epsilon=array2D ( SII+1,SJJ+1) ;
//21. dimensionless , relative dielectric constant
double **mu=array2D ( SII+1,SJJ+1) ;
//22. dimensionless , relative permeability
double *aax=array1D ( SII+1) ; double *bbx=array1D ( SII+1) ;
131
double *ccx=array1D ( SII+1) ;//23 CPML coefficients
double *aay=array1D ( SJJ+1) ; double *bby=array1D ( SJJ+1) ;
double *ccy=array1D ( SJJ+1) ;//24 CPML coefficients
// double *a_Hz=array1D ( SII+1) ; double *b_Hz=array1D ( SII+1) ;
// double *c_Hz=array1D ( SII+1) ;// CPML coefficients
double *kkx=array1D ( SII+1) ; double *kky=array1D ( SJJ+1) ;
// double *k_Hz=array1D ( SII+1) ;//25 CPML coefficients
double *sigmax=array1D ( SII+1) ; / / 2 6 . conductivity component , CPML terms
double *sigmay=array1D ( SJJ+1) ; / / 2 7 . conductivity component
// double ** sigmastarx=array2D ( SII+1,SJJ+1) ;
//27 . 1 . magnetic loss component ( fictional )
double *** psi_dyHz=array3D (3 , SII+1,SJJ+1) ; / / 2 8 . CPML
double *** psi_dxHz=array3D (3 , SII+1,SJJ+1) ; / / 2 9 . CPML
double *** psi_dxEy=array3D (3 , SII+1,SJJ+1) ; / / 3 0 . CPML
double *** psi_dyEx=array3D (3 , SII+1,SJJ+1) ; / / 3 1 . CPML
double *dx=array1D ( SII+1) ;
//32. units : nm , grids s i z e , 1/ Nlambday of wavelength 500nm ,
// typical optical value (500 nm ) ˜ 25nm
double *dy=array1D ( SJJ+1) ;//33
double *hxi=array1D ( SII+1) ;//340 dual edge
double *hyj=array1D ( SJJ+1) ;//35
double *xx=array1D ( SII+1) ;//36 abs x−ax i s value
double *yy=array1D ( SJJ+1) ;//37
double * HsourceAmp=array1D ( SII+1) ; / / 3 8 . output terms
double *EHphase=array1D ( SII+1) ; / / 3 9 . shift phase to match E
double ** HzAmp=array2D ( SII+1,SJJ+1) ;
//40. Average intensity of Hz , 2 periods average after being in steady state
double ** ExAmp=array2D ( SII+1,SJJ+1) ; / / 4 1 .
double ** EyAmp=array2D ( SII+1,SJJ+1) ; / / 4 2 .
132
double *** ET2=array3D (3 , SII+1,SJJ+1) ; //43 . 4 | E |ˆ2
double ** DetectorsR4=array2D (6* ( int ) nT+2,SJJ+1) ;
// Detectors in r e a l time one Rightside 4/4 of SII
double ** DetectorsT=array2D (6* ( int ) nT+2,SII+1) ;// Detectors on Top
double ** DetectorsB=array2D (6* ( int ) nT+2,SII+1) ;// Detectors on Bottom
double ** DetectorsR1=array2D (6* ( int ) nT+2,SJJ+1) ;
// Detectors in r e a l time one Rightside 1/4 of SII position
double ** DetectorsR2=array2D (6* ( int ) nT+2,SJJ+1) ;
// Detectors in r e a l time one Rightside 2/4 of SII position
double ** DetectorsR3=array2D (6* ( int ) nT+2,SJJ+1) ;
// Detectors in r e a l time one Rightside 3/4 of SII position
FullCubicFormula cubicRoot ; // to solve f o r E from D analytically
FILE *fp_Hz ; /* declare a File Pointer f o r Hz */
FILE *fp_Ey ; /* declare a File Pointer f o r Ey */
FILE *fp_Ex ;
FILE* fp_ADE ;//1−D auxiliary grids ;
FILE *fp_info ; // additional i n f o file
FILE* fp_Spectrum ;
short file_importflag=1;//1 f o r import data to continue
short file_exportflag=0;//1 f o r ready to export a l l final fields data
short file_fullflag=0;//1 f o r additional files output
double Hup , Hdown , Hleft , Hright ; // f o r r e a l values : up , down , right , left
double Eup , Edown , Eleft , Eright ;
// substitution f o r Hz and E ' s in Yee formulae ,
// especially f o r TFSF_n calculation
int distance_n =0;
// distance in terms of cell number , from a point to bottom left corner of TF
133
/*** Startline of main program body−−−−−−−−−−−−−−−−−−−−−−*/
/***** Startline of ALL cell value initialization−−−−−−−*/
/******* Startline of nonuniform cell generator−−−−−−−−*/
i f ( ( gradeRatey >=1.4) | | ( gradeRatey <=0.6) ){
cout<<gradeRatey<<endl ;
cout<<
”Varying cell step too big ! Verify parameters f o r \
non−uniform grids !\ n\n ” ;
cout<<”Press 'q ' to quit , or others to continue .\ n\n ” ;
char gradequit ;
cin>>gradequit ;
i f ( gradequit== 'q ' )
exit (0) ;
}
f o r ( ii=0;ii<=SII ; ii++){//dx , hxi initialization to uniform g r i d value
dx [ ii ]=dxmax ; // dx , dy are dividing Hz field
}
f o r ( jj=0;jj<=SJJ ; jj++){//dy , hyj initialization to uniform g r i d value
dy [ jj ]=dymax ;
}
/* i f (0 ) {// nonuniform g r i d mesh , x
f o r ( ii=gradeLeft ; ii<gradeMinLeft ; ii++){// grade down
dx [ ii ]=dxmax *pow ( gradeRatex , ii−gradeLeft ) ;
}
f o r ( ii=gradeMinLeft ; ii<gradeMinRight ; ii++){//min basin
dx [ ii ]=dxmin ;
}
f o r ( ii=gradeMinRight ; ii<gradeRight ; ii++){// grade up
dx [ ii ]=dxmax *pow ( gradeRatex , gradeRight−ii ) ;
}
134
}// i f , x
*/ int levelYcenter=centerJ ; // basin valley center
int levely0=levelYcenter−mingrade_n−varygrade_n ; // varying low limit
int levely1=levelYcenter−mingrade_n ; // basin low edge
int levely2=levelYcenter+mingrade_n ; // basin high edge
int levely3=levelYcenter+mingrade_n+varygrade_n ; // varying high limit
i f (1 ) {// nonuniform grids , y
f o r ( jj=levely0 ; jj<levely1 ; jj++){// low , dy decreasing
dy [ jj ]=dymax *pow ( gradeRatey , levely0−jj ) ;
// cout<<jj<<”\t”<<dy [ jj]<<endl ;
}
f o r ( jj=levely1 ; jj<levely2 ; jj++){// middle , , dy const
dy [ jj ]=dymin ;
// cout<<jj<<”\t”<<dy [ jj]<<endl ;
}
f o r ( jj=levely2 ; jj<levely3 ; jj++){// high , dy increasing
dy [ jj ]=dymax *pow ( gradeRatey , jj−levely3 ) ;
// cout<<jj<<”\t”<<dy [ jj]<<endl ;
}
}// i f y
f o r ( ii=1;ii<=SII ; ii++){//dx , hxi initialization to uniform g r i d value
hxi [ ii ]=( dx [ ii ]+dx [ ii−1]) /2;// hxi , hyj are dividing Ex , Ey fields
xx [ ii ]=xx [ ii−1]+dx [ ii ] ;
}
hxi [ 0 ]= dx [ 0 ] ;
f o r ( jj=1;jj<=SJJ ; jj++){//dy , hyj initialization to uniform g r i d value
hyj [ jj ]=( dy [ jj ]+dy [ jj−1]) /2 ;
yy [ jj ]=yy [ jj−1]+dy [ jj ] ;
}
hyj [ 0 ]= dy [ 0 ] ;
135
double dxTFl=dx [ TFl ] ; / / f o r Hz correction
double hxiTFl_1=hxi [ TFl −1] ;// f o r Dy correction
/*−−end of nonuniform cell generator ***********************/
/***** start of CPML parameters−−−−−−−−−−−−−−−−−−−−−−−*/
double sigma_max =0.0005;//0.08 −0.0003
double kkx_max =1.50;//s=1/(k_w+sigma_w /( a_w+j omega epsilon_0 ) ) , w=x , y , z
double aax_max =0.001;//keep this term zero f o r most cases !
f o r ( ii=0;ii<=SII ; ii++){
i f ( ii<=SFl ) {// CPML left side
double powerleft=pow ( ( double ) ( SFl−ii ) /CPML_n , 2 ) ;
kkx [ ii ]=1+(kkx_max −1)*powerleft ;
sigmax [ ii ]=powerleft * sigma_max ;
aax [ ii ]=aax_max *pow ( ( ( double ) ii ) /CPML_n , 2 ) ;
bbx [ ii ]=exp(−dt *( sigmax [ ii ] / epsilon_0 /kkx [ ii ]+aax [ ii ] / epsilon_0 ) ) ;
ccx [ ii ]=( bbx [ ii ]−1) /( kkx [ ii ]+kkx [ ii ]* kkx [ ii ]* aax [ ii ] / sigmax [ ii ] ) ;
i f ( sigmax [ ii ]==0)
ccx [ ii ]=0;
}
e l s e i f ( ii>=SFr ) {// CPML right side
double powerright=pow ( ( double ) (ii−SFr ) /CPML_n , 2 ) ;
sigmax [ ii ]=powerright *sigma_max ;
kkx [ ii ]=1+(kkx_max −1)*powerright ;
aax [ ii ]=aax_max *pow ( ( double ) ( SII−ii ) /CPML_n , 2 ) ;
bbx [ ii ]=exp(−dt *( sigmax [ ii ] / epsilon_0 /kkx [ ii ]+aax [ ii ] / epsilon_0 ) ) ;
ccx [ ii ]=( bbx [ ii ]−1) /( kkx [ ii ]+kkx [ ii ]* kkx [ ii ]* aax [ ii ] / sigmax [ ii ] ) ;
i f ( sigmax [ ii ]==0)
ccx [ ii ]=0;
}
e l s e {//0
sigmax [ ii ]=0;
kkx [ ii ]=1;
136
aax [ ii ]=0;
}
}// f o r ii
double kkxTFl=kkx [ TFl ] ; / / f o r Hz , E correction
f o r ( jj=0;jj<=SJJ ; jj++){
i f ( jj<=SFb ) {// CPML bottom
double powerbottom=pow ( ( double ) ( SFb−jj ) /CPML_n , 2 ) ;
sigmay [ jj ]=powerbottom *sigma_max ;
kky [ jj ]=1+(kkx_max −1)*powerbottom ;
aay [ jj ]=aax_max *pow ( ( ( double ) ( SFb−1)−jj ) /CPML_n , 2 ) ;
bby [ jj ]=exp(−dt *( sigmay [ jj ] / epsilon_0 /kky [ jj ]+aay [ jj ] / epsilon_0 ) ) ;
ccy [ jj ]=( bby [ jj ]−1) /( kky [ jj ]+kky [ jj ]* kky [ jj ]* aay [ jj ] / sigmay [ jj ] ) ;
i f ( sigmay [ jj ]==0)
ccy [ jj ]=0;
}
e l s e i f ( jj>=SFt ) {// CPML top
double powertop=pow ( ( double ) ( jj−SFt ) /CPML_n , 2 ) ;
sigmay [ jj ]=powertop *sigma_max ;
kky [ jj ]=1+(kkx_max −1)*powertop ;
aay [ jj ]=aax_max *pow ( ( double ) ( SJJ−jj ) /CPML_n , 2 ) ;
bby [ jj ]=exp(−dt *( sigmay [ jj ] / epsilon_0 /kky [ jj ]+aay [ jj ] / epsilon_0 ) ) ;
ccy [ jj ]=( bby [ jj ]−1) /( kky [ jj ]+kky [ jj ]* kky [ jj ]* aay [ jj ] / sigmay [ jj ] ) ;
i f ( sigmay [ jj ]==0)
ccy [ jj ]=0;
}
e l s e {
sigmay [ jj ]=0;
kky [ jj ]=1;
aay [ jj ]=0;
}
}// f o r jj
137
double chi3 [ SII+1]={0};
f o r ( ii=TFl ; ii<=SII ; ii++){
chi3 [ ii ]=chi3M ;
//(1− exp(−(ii−TFl ) *( ii−TFl ) /1.0/ Nlambdax /Nlambdax ) ) *chi3M ;
// chi3 ramp up
}
f o r ( ii=0;ii<=SII ; ii++){
f o r ( jj=0;jj<=SJJ ; jj++){
epsilon [ ii ] [ jj ]=1 .0 ; // medium initialization , relative values
mu [ ii ] [ jj ]=1 . 0 ;
}// f o r ( jj )
}// f o r ( ii )
/*−−−−−end of CPML parameters ************************/
/*−−−−−Endline of ALL cell initialization ***********/
/****** Startline of FDTD major loops−−−−−−−−−−*/
double smoothupHc =0.001;
// unit : V/m artificial factor to tune up H field slowly
double smoothupEc=smoothupHc ;
// together with chi (3) m2/v2 , the intensity unit is W/m2
double smoothupH =0;// ini .
double smoothupE =0;// ini .
double decayfacy =1;// ini . along y−axi s , decay term f o r SP
double phase_kx=atan2 ( kx_I , kx_R ) ;
// phase difference due to wavevector Qx between H and E fields
double abs_kx=sqr t ( kx_I*kx_I+kx_R*kx_R ) ;
// absolute value of wavevector Qx
short rampup=4;// total time steps to ramp up the fields
double *chi3ii ; // nonlinear susceptibility
double *Exnp1 ,* Hznp1 ,* Eynp1 ; // temp variables to speed up loop calculation
double *Exnp2 ,* Hznp2 ,* Eynp2 ; // temp variables to speed up loop calculation
138
double *Dxnp2 ,* Dynp2 ; // temp variables to speed up loop calculation
double *Px_Drudenp2 ,* Py_Drudenp2 ;
double *ET2np1 ,* ET2np2 ;
double *psi_dyHznp1 ,* psi_dxHznp1 ;
double *psi_dyHznp2 ,* psi_dxHznp2 ;
double *hyjjj ,* hxiii ,* dxii ,* dyjj ;
double *Px_Kerrnp1 ,* Py_Kerrnp1 ;
double *psi_dxEynp0 ,* psi_dyExnp0 ;
double yymetalJ=yy [ metalJ ] ;
//*************************************************
f o r ( nn=1;nn<=Ntotal ; nn++){// time loop
np0=(nn−1)%3 ; np1=nn%3;np2=(nn+1)%3;// pa r i t y convert ion , 3−based f o r ADE
smoothupH=smoothupHc *(1−exp(−(double ) nn/rampup /nT ) ) ;
// f o r 1W laser on 2mmx2mm spot , field strength ˜1000 V/m
smoothupE=smoothupEc *(1−exp(−(nn+0.5)/rampup /nT ) ) ;
//Boyd , Nonlinear Optics , Page3 E=5.14 x10 ˆ11 V/m
f o r ( jj=0;jj<=metalJ ; jj++){// SPP source , in Metal
ydis=yy [ jj ]−yymetalJ ; // ydis<0
decayfacy=exp ( ky2_R*ydis ) ;
//Hz and Ey have the same y coordinates position
Hzin [ nn%2 ] [ j j ]=smoothupH* s i n ( omega 0*nn*dt ) *decayfacy ;
//Hz , at step n , position i//
Eyin [ nn%2 ] [ j j ]=smoothupE* s i n ( omega 0 *( nn+0.5)*dt+kx R*dxTFl*0.50−
phase_kx ) * decayfacy *abs_kx /omega_0 /epsilon_0 / epsilon_inf ;
//Dy , at step n+1/2 , position i−1/2
}
f o r ( jj=metalJ+1;jj<=SJJ ; jj++){// SPP source , in Air
ydis=yy [ jj ]−yymetalJ ; //
decayfacy=exp(−ky0_R*ydis ) ;
Hzin [ nn%2 ] [ j j ]=smoothupH* s i n ( omega 0*nn*dt ) *decayfacy ;
139
//Hz , at step n , position i
Eyin [ nn%2 ] [ j j ]=smoothupE* s i n ( omega 0 *( nn+0.5)*dt+kx R*dxTFl*0.50−
phase_kx ) * decayfacy *abs_kx /omega_0 /epsilon_0 ;
//Dy , at step n+1/2 , position i−1/2
}
f o r ( ii=0;ii<=SII−1;ii++){// Terms w/o E [ n+1]
f o r ( jj=0;jj<=SJJ−1;jj++){
Hup=Hz [ np2 ] [ ii ] [ jj+1] ;//H , named relative to E
Hdown=Hz [ np2 ] [ ii ] [ jj ] ;
Hright=Hz [ np2 ] [ ii+1] [ jj ] ; / / H , named relative to E
Hleft=Hz [ np2 ] [ ii ] [ jj ] ;
chi3ii=&chi3 [ ii ] ;
ET2np1=&ET2 [ np1 ] [ ii ] [ jj ] ;
Exnp1=&Ex [ np1 ] [ ii ] [ jj ] ;
Eynp1=&Ey [ np1 ] [ ii ] [ jj ] ;
Px_Drudenp2=&Px_Drude [ np2 ] [ ii ] [ jj ] ;
Py_Drudenp2=&Py_Drude [ np2 ] [ ii ] [ jj ] ;
Px_Kerrnp1=&Px_Kerr [ np1 ] [ ii ] [ jj ] ;
Py_Kerrnp1=&Py_Kerr [ np1 ] [ ii ] [ jj ] ;
/*−−−Polarization calc starts−−−−−−−−−−−−−−*/
i f ( ( jj<=metalJ ) ) {// in metal
//−−−−−−terms include E [ n]−−−−−−−−−−−−−−
* Px_Drudenp2=a_drude *Px_Drude [ np1 ] [ ii ] [ jj ]+b_drude *
Px_Drude [ np0 ] [ ii ] [ jj ]+c_drude *(* Exnp1 ) ;
// c_drude includes Epsilon_0 already
* Py_Drudenp2=a_drude *Py_Drude [ np1 ] [ ii ] [ jj ]+b_drude *
Py_Drude [ np0 ] [ ii ] [ jj ]+c_drude *(* Eynp1 ) ;
// Px_Lorentz [ np2 ] [ ii ] [ jj ]=0;
// a_lorentz *Px_Lorentz [ np1 ] [ ii ] [ jj ]+b_lorentz *
// Px_Lorentz [ np0 ] [ ii ] [ jj ]+c_lorentz *Ex [ np1 ] [ ii ] [ jj ] ;
// Py_Lorentz [ np2 ] [ ii ] [ jj ]=0;
140
// a_lorentz *Py_Lorentz [ np1 ] [ ii ] [ jj ]+b_lorentz
//* Py_Lorentz [ np0 ] [ ii ] [ jj ]+c_lorentz *Ey [ np1 ] [ ii ] [ jj ] ;
// Sp [ np2 ] [ ii ] [ jj ]=0;// a_raman *Sp [ np1 ] [ ii ] [ jj ]+b_raman *
//Sp [ np0 ] [ ii ] [ jj ]+c_raman *( Ex [ np1 ] [ ii ] [ jj ]* Ex [ np1 ] [ ii ] [ jj ]+
// Ey [ np1 ] [ ii ] [ jj ]* Ey [ np1 ] [ ii ] [ jj ] ) ;
//−−−−−−−−−−−terms include E [ n+1]−−−−−−−−−−−−−−−−−−−−−−
// Px_Debye [ np1 ] [ ii ] [ jj ]=0;
// a_debye * Px_Debye [ np0 ] [ ii ] [ jj ]+b_debye *( Ex [ np1 ] [ ii ] [ jj ]+
//Ex [ np0 ] [ ii ] [ jj ] ) ;
// Py_Debye [ np1 ] [ ii ] [ jj ]=0;
// a_debye * Py_Debye [ np0 ] [ ii ] [ jj ]+b_debye *( Ey [ np1 ] [ ii ] [ jj ]+
//Ey [ np0 ] [ ii ] [ jj ] ) ;
* Px_Kerrnp1=alpha*epsilon_0 *(* chi3ii ) *(* ET2np1 ) *(* Exnp1 ) ;
// epsilon_0 included
* Py_Kerrnp1=alpha*epsilon_0 *(* chi3ii ) *(* ET2np1 ) *(* Eynp1 ) ;
// Px_Raman [ np1 ] [ ii ] [ jj ]=epsilon_0 *Ex [ np1 ] [ ii ] [ jj ]*
//Sp [ np1 ] [ ii ] [ jj ] ;
// Py_Raman [ np1 ] [ ii ] [ jj ]=epsilon_0 *Ey [ np1 ] [ ii ] [ jj ]*
//Sp [ np1 ] [ ii ] [ jj ] ;
}// metal
e l s e {// in air
//−−−−−−−−−−terms include E [ n ]
// * Px_Drudenp2 =0;
// * Py_Drudenp2 =0;
// Px_Lorentz [ np2 ] [ ii ] [ jj ]=0;
// Py_Lorentz [ np2 ] [ ii ] [ jj ]=0;
// Sp [ np2 ] [ ii ] [ jj ]=0;
//−−−−−−−−−−−terms include E [ n+1]
// Px_Debye [ np1 ] [ ii ] [ jj ]=0;
// Py_Debye [ np1 ] [ ii ] [ jj ]=0;
// * Px_Kerrnp1 =0;
141
// * Py_Kerrnp1 =0;
// Px_Raman [ np1 ] [ ii ] [ jj ]=0;
// Py_Raman [ np1 ] [ ii ] [ jj ]=0;
}// air
}
}
/*−−−−−−−Polarization calc ends−−−−−−−−−−−−−−−−−−−−−*/
/*−−−−−−−−−Displacement field starts−−−−−−−−*/
f o r ( ii=0;ii<=SII ; ii++){// Dx loop−−−−−−−−−−−−−−−−
f o r ( jj=0;jj<=SJJ−1;jj++){
Hup=Hz [ np2 ] [ ii ] [ jj+1] ;//H , named relative to E
Hdown=Hz [ np2 ] [ ii ] [ jj ] ;
psi_dyHznp1=&psi_dyHz [ np1 ] [ ii ] [ jj ] ;
psi_dyHznp2=&psi_dyHz [ np2 ] [ ii ] [ jj ] ;
hyjjj=&hyj [ jj ] ;
*psi_dyHznp2=bby [ jj ]* (* psi_dyHznp1 )+ccy [ jj ] * ( Hup−Hdown ) /(* hyjjj ) ;
// CPML term , 0 unless in CPML layer
Dx [ np2 ] [ ii ] [ jj ]=Dx [ np1 ] [ ii ] [ jj ]+dt*
( ( Hup−Hdown ) /(* hyjjj ) /kky [ jj ]+(* psi_dyHznp2 ) ) ;
}// f o r ( jj )
}//E loop ( ii )
f o r ( ii=0;ii<=SII−1;ii++){// Dy loop−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
f o r ( jj=0;jj<=SJJ ; jj++){
Hright=Hz [ np2 ] [ ii+1] [ jj ] ; / / H , named relative to E
Hleft=Hz [ np2 ] [ ii ] [ jj ] ;
psi_dxHznp2=&psi_dxHz [ np2 ] [ ii ] [ jj ] ;
psi_dxHznp1=&psi_dxHz [ np1 ] [ ii ] [ jj ] ;
hxiii=&hxi [ ii ] ;
142
*psi_dxHznp2=bbx [ ii ]* (* psi_dxHznp1 )+
ccx [ ii ] * ( Hright−Hleft ) /(* hxiii ) ; // CPML term
Dy [ np2 ] [ ii ] [ jj ]=Dy [ np1 ] [ ii ] [ jj ]−dt*
( ( Hright−Hleft ) /(* hxiii ) /kkx [ ii ]+(* psi_dxHznp2 ) ) ;
}// f o r ( jj )
}//E loop ( ii )
f o r ( jj=SFb ; jj<=SFt ; jj++){// Dy TFSF source correction , Hz : TF−>SF
Hright=Hz [ np2 ] [ TFl ] [ jj ]−Hzin [ nn%2 ] [ j j ] ; / /H, named r e l a t i v e to E
Hleft=Hz [ np2 ] [ TFl −1] [ jj ] ;
Dy [ np2 ] [ TFl −1] [ jj ]=Dy [ np1 ] [ TFl −1] [ jj ]+dt *( Hleft−Hright ) / hxiTFl_1 ;
// Dy [ np2 ] [ TFl −1] [ jj ]=( exp(−tau *0.5* dt ) *Dy [ np1 ] [ TFl −1] [ jj ]+
//dt *( Hleft−Hright ) /hxi [ TFl−1]) /( exp ( tau *0.5* dt ) ) ;
}
f o r ( ii=0;ii<=SII ; ii++){
//Ex Ey loop−−−−−−−−t=n+1/2−−−−−−−−−−−−−−−−−−−−−−−−−−−−
f o r ( jj=0;jj<=SJJ ; jj++){
Exnp2=&Ex [ np2 ] [ ii ] [ jj ] ; / / to be calc
Eynp2=&Ey [ np2 ] [ ii ] [ jj ] ; / / to be calc
Dxnp2=&Dx [ np2 ] [ ii ] [ jj ] ; / / f i x , from H
Dynp2=&Dy [ np2 ] [ ii ] [ jj ] ; / / f i x , from H
i f ( Metalflag ( ii , jj ) ){
ET2np2=&ET2 [ np2 ] [ ii ] [ jj ] ;
// ET2np1=&ET2 [ np1 ] [ ii ] [ jj ] ;
Px_Drudenp2=&Px_Drude [ np2 ] [ ii ] [ jj ] ;
Py_Drudenp2=&Py_Drude [ np2 ] [ ii ] [ jj ] ;
chi3ii=&chi3 [ ii ] ;
double cubic_a0 , cubic_a , cubic_b , cubic_c , AA ; // AA=E*E
cubic_a0=square (* chi3ii ) ;
cubic_a=2*epsilon_inf *(* chi3ii ) ;
cubic_b=epsilon_inf *epsilon_inf ;
143
cubic_c=−(square (* Dxnp2−*Px_Drudenp2 )+square (* Dynp2−
*Py_Drudenp2 ) ) /( epsilon_0 *epsilon_0 ) ;
cubicRoot . setFormula ( cubic_a0 , cubic_a , cubic_b , cubic_c ) ;
AA=r ea l ( cubicRoot . cubicRoot1 ( ) ) ;
*ET2np2=AA ;
// cout<<”a , b , c\t”<<cubic_a<<cubic_b<<cubic_c<<endl ;
// cout<<AA<<endl ;
*Exnp2=((* Dxnp2 )−(*Px_Drudenp2 ) ) /epsilon_0 /
( epsilon_inf+(*chi3ii ) *AA ) ;
*Eynp2=((* Dynp2 )−(*Py_Drudenp2 ) ) /epsilon_0 /
( epsilon_inf+(*chi3ii ) *AA ) ;
// *ET2np2=*ET2np1 ;
// f o r ( int mm=1;mm<=nonlinearM ; mm++){
//Ex nonlinear iteration , f o r P (w=w+w−w ) only
// double Dynp2compare ;
// Dynp2compare=*Dynp2 ;
//// f o r ( mm=0;fabs ( Dynp2compare−*Dynp2 ) >=0.0000; mm++){
// *Exnp2=(*Dxnp2/*−a_debye *Px_Debye [ np1 ] [ ii ] [ jj ]−
// b_debye *Ex [ np1 ] [ ii ] [ jj ]−Px_Lorentz [ np2 ] [ ii ] [ jj ]*/−
//(* Px_Drudenp2 ) ) /( epsilon_0 *epsilon_inf /*+b_debye*/+
//alpha *epsilon_0 *(* chi3ii ) *(* ET2np2 )/*+
// epsilon_0 *Sp [ np2 ] [ ii ] [ jj ]*/ ) ;
// *Eynp2=(*Dynp2/*−a_debye *Py_Debye [ np1 ] [ ii ] [ jj ]−
// b_debye *Ey [ np1 ] [ ii ] [ jj ]−Py_Lorentz [ np2 ] [ ii ] [ jj ]*/−
//(* Py_Drudenp2 ) ) /( epsilon_0 *epsilon_inf /*+
// b_debye*/+alpha*epsilon_0 *(* chi3ii ) *(* ET2np2 )/*+
// epsilon_0 *Sp [ np2 ] [ ii ] [ jj ]*/ ) ;
// *ET2np2=(*Exnp2 ) *(* Exnp2 )+(*Eynp2 ) *(* Eynp2 ) ;
// Dynp2compare=*Eynp2 *( epsilon_0 * epsilon_inf+
//alpha *epsilon_0 *(* chi3ii ) *(* ET2np2 ) )+
//* Py_Drudenp2 ;
//compared , back calc value
144
// i f ( ( nn>500)&&(jj==metalJ )&&(ii==TFl+10) ){
// cout<<”nn , ii , jj==”<<nn<<”\t”<<ii<<”\t”<<jj<<endl ;
// // system (” pause ”) ;
// cout<<”Ex , Ey”<<*Exnp2<<”\t”<<*Eynp2<<endl ;
// cout<<”Dynp2compare==”<<Dynp2compare<<endl ;
// cout<<”*Dynp2==”<<*Dynp2<<endl ;
// cout<<”Differ==”<<Dynp2compare−*Dynp2<<”\tOR \t”
//<<(Dynp2compare /(* Dynp2 )−1)*100<<”%”<<endl<<endl ;
// }// i f nn
// }// f o r ( mm iteration )
}// i f metalflag
e l s e i f ( Airflag (ii , jj ) ) {// in air
*Exnp2=(*Dxnp2 ) /epsilon_0 ;
*Eynp2=(*Dynp2 ) /epsilon_0 ;
}// e l s e i f Airflag
}// f o r ( jj )
}//E loop ( ii )
f o r ( ii=1;ii<=SII−1;ii++){// Hz loop−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
f o r ( jj=1;jj<=SJJ−1;jj++){
Eright=Ey [ np2 ] [ ii ] [ jj ] ;
Eleft=Ey [ np2 ] [ ii−1] [ jj ] ;
Eup=Ex [ np2 ] [ ii ] [ jj ] ;
Edown=Ex [ np2 ] [ ii ] [ jj−1] ;
psi_dxEynp0=&psi_dxEy [ np0 ] [ ii ] [ jj ] ;
psi_dyExnp0=&psi_dyEx [ np0 ] [ ii ] [ jj ] ;
dxii=&dx [ ii ] ;
dyjj=&dy [ jj ] ;
*psi_dxEynp0=bbx [ ii ]* psi_dxEy [ np2 ] [ ii ] [ jj ]+
145
ccx [ ii ] * ( Eright−Eleft ) /(* dxii ) ; // CPML term
*psi_dyExnp0=bby [ jj ]* psi_dyEx [ np2 ] [ ii ] [ jj ]+
ccy [ jj ] * ( Eup−Edown ) /(* dyjj ) ; // CPML term
Hz [ np0 ] [ ii ] [ jj ]=Hz [ np2 ] [ ii ] [ jj ]−
dt/mu_0 *((1/ kkx [ ii ] * ( Eright−Eleft ) /(* dxii )−
1/kky [ jj ] * ( Eup−Edown ) /(* dyjj )+*psi_dxEynp0−*psi_dyExnp0 ) ) ;
}// f o r ( jj )
}// Hz loop ( ii )
i f (1 ) {// Hz correction
f o r ( jj=SFb ; jj<=SFt ; jj++){// Hzx TFSF correction , Ey : SF−>TF
Eright=Ey [ np2 ] [ TFl ] [ jj ] ;
Eleft=Ey [ np2 ] [ TFl −1] [ jj ]+Eyin [ nn%2 ] [ j j ] ;
Eup=Ex [ np2 ] [ TFl ] [ jj ] ;
Edown=Ex [ np2 ] [ TFl ] [ jj−1] ;
Hz [ np0 ] [ TFl ] [ jj ]=Hz [ np2 ] [ TFl ] [ jj ]−
dt/mu_0 *((1/ kkxTFl *( Eright−Eleft ) /dxTFl−
1/kky [ jj ] * ( Eup−Edown ) /dy [ jj ]+psi_dxEy [ np0 ] [ TFl ] [ jj ]−
psi_dyEx [ np0 ] [ TFl ] [ jj ] ) ) ;
}
}
double *ExAmpii ,* EyAmpii ,* HzAmpii ;
i f ( nn>Ntotal−2*nT ){
// amplitude search , 2 period of time whi le being in steady state
f o r ( ii=0;ii<=SII ; ii++){
f o r ( jj=0;jj<=SJJ ; jj++){
ExAmpii=&ExAmp [ ii ] [ jj ] ;
EyAmpii=&EyAmp [ ii ] [ jj ] ;
HzAmpii=&HzAmp [ ii ] [ jj ] ;
// HsourceAmp [ jj ]=max( HsourceAmp [ jj ] , abs ( Hzin [ ( nn+1)%2 ] [ j j ] ) ) ;
146
*ExAmpii=max( fabs ( Ex [ np2 ] [ ii ] [ jj ] ) ,* ExAmpii ) ;
*EyAmpii=max( fabs ( Ey [ np2 ] [ ii ] [ jj ] ) ,* EyAmpii ) ;
*HzAmpii=max( fabs ( Hz [ np0 ] [ ii ] [ jj ] ) ,* HzAmpii ) ;
}// jj
}// ii
}// amplitude search loop ends
i f ( nn>Ntotal −6*(int ) nT ) {// spectrum data detector
f o r ( jj=0;jj<=SJJ ; jj++){
DetectorsR4 [ nn−Ntotal+6*(int ) nT ] [ jj ]=Hz [ np0 ] [ TFr −1] [ jj ] ;
DetectorsR1 [ nn−Ntotal+6*(int ) nT ] [ jj ]=Hz [ np0 ] [ SII / 4 ] [ jj ] ;
DetectorsR2 [ nn−Ntotal+6*(int ) nT ] [ jj ]=Hz [ np0 ] [ SII / 2 ] [ jj ] ;
DetectorsR3 [ nn−Ntotal+6*(int ) nT ] [ jj ]=Hz [ np0 ] [ SII * 3 / 4 ] [ jj ] ;
}
f o r ( ii=0;ii<=SII ; ii++){
DetectorsT [ nn−Ntotal+6*(int ) nT ] [ ii ]=Hz [ np0 ] [ ii ] [ SFt −1] ;
DetectorsB [ nn−Ntotal+6*(int ) nT ] [ ii ]=Hz [ np0 ] [ ii ] [ SFb+1] ;
}
}
/*−−−−−Endline of FDTD major loops *************************/
/******* Startline of . dat file writing−−−−−−−−−−−−−−−*/
static clock_t lasttimeclock=clock ( ) ;
i f ( nn%10==0) {
cout<<”Time Step = ”<<nn<<endl ; // step screen display monitor
i f ( nn%50==0){
timer_yan ( Ntotal , nn , 5 0 , lasttimeclock ) ;
lasttimeclock=clock ( ) ;
}
}
FileFlag=((nn%f i l e s t e p==0)&&(nn!=0) ) ;
147
char namenumber [ 1 8 ] ; / / extension of output files
i f ( FileFlag ) {// Hz auto file writing
_itoa ( nn , namenumber , 1 0 ) ;
// convert time step into file name string , 10 based
strcat ( namenumber , ” . dat ”) ;
char file_Hz [ 18 ]=” Hz ” ;
strcat ( file_Hz , namenumber ) ; // define output ( data ) Hz
fp_Hz = fopen ( file_Hz , ”w ”) ;
// f p r i n t f ( fp_Hz , ” x−ax i s \ty−ax i s \tHz\n ”) ;
// f p r i n t f ( fp_Hz , ” nm\tnm\tA . U .\ n ”) ;
//Hz , column out put starts
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii+=2){
f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj+=2){
f p r i n t f ( fp_Hz , ”%f \ t%f \ t%f \n” , xx [ i i ] , yy [ j j ] , Hz [ np2 ] [ i i ] [ j j ] ) ;
// column output format
}// f o r jj
}// f o r ii//Hz column output ends
//Hz output , matrix output format starts
/* f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj++){//
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii++)
{
f p r i n t f ( fp_Hz , ”%f \ t ” ,Hz [ np2 ] [ i i ] [ j j ] ) ;
}// f o r ii
f p r i n t f ( fp_Hz , ”\ n ”) ;
}// f o r jj// matrix output end
*/ printf (”\ t\t\tFile written :%s \n” , f i l e H z ) ;
//step indicator on screen
f c l o s e ( fp_Hz ) ;
}// file Hz writing ends ( i f )
148
i f ( FileFlag ) {// Ey auto file writing
itoa (nn , namenumber , 1 0 ) ;
// convert time step into file name string , 10 based
strcat ( namenumber , ” . dat ”) ;
char file_Ey [ 18 ]=” Ey ” ;
strcat ( file_Ey , namenumber ) ; // define output ( data ) Ey
fp_Ey = fopen ( file_Ey , ”w ”) ;
// f p r i n t f ( fp_Ey , ” x−ax i s \ty−ax i s \tEy\n ”) ;
// f p r i n t f ( fp_Ey , ” nm\tnm\tA . U .\ n ”) ;
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii+=2){
f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj+=2){
f p r i n t f ( fp_Ey , ”%l f \ t%l f \ t%l f \n” , xx [ i i ] , yy [ j j ] , Ey [ np2 ] [ i i ] [ j j ] ) ;
// column output format
}
}
// f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj++){// Ey output
// f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii++)
// {
// f p r i n t f ( fp_Ey , ”%f \ t ” ,Ey [ np2 ] [ i i ] [ j j ] ) ;
// }
// f p r i n t f ( fp_Ey , ”\ n ”) ;
//}
f c l o s e ( fp_Ey ) ;
}// file Ey writing ends ( i f )
i f ( FileFlag ) {// Ex auto file writing
itoa (nn , namenumber , 1 0 ) ;
149
// convert time step into file name string , 10 based
strcat ( namenumber , ” . dat ”) ;
char file_Ex [ 18 ]=” Ex ” ;
strcat ( file_Ex , namenumber ) ; // define output ( data ) Ex
fp_Ex = fopen ( file_Ex , ”w ”) ;
// f p r i n t f ( fp_Ex , ” x−ax i s \ty−ax i s \tEx\n ”) ;
// f p r i n t f ( fp_Ex , ” nm\tnm\tA . U .\ n ”) ;
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii+=2){
f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj+=2){
f p r i n t f ( fp_Ex , ”%l f \ t%l f \ t%l f \n” , xx [ i i ] , yy [ j j ] , Ex [ np2 ] [ i i ] [ j j ] ) ;
// column output format
}
}
// f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj++){// Ex output
// f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii++)
// {
// f p r i n t f ( fp_Ex , ”%f \ t ” ,Ex [ np2 ] [ i i ] [ j j ] ) ;
// }
// f p r i n t f ( fp_Ex , ”\ n ”) ;
//}
f c l o s e ( fp_Ex ) ;
}// file Ex writing ends ( i f )
i f (0 ) {// output dx stucture
fp_Hz=fopen (” dxminmax . dat ” ,” w ”) ;
f o r ( ii=0;ii<=SII−1;ii++){
f p r i n t f ( fp_Hz , ”%d\ t%l f \n” , i i , dx [ i i ] ) ;
}
150
f c l o s e ( fp_Hz ) ;
}
i f ( FileFlag&&0) {// ADE file writing
itoa (nn , namenumber , 1 0 ) ;
// convert time step into file name string , 10 based
strcat ( namenumber , ” . dat ”) ;
char file_ADE [ 18 ]=” ADE_Hz ” ;
strcat ( file_ADE , namenumber ) ;
fp_ADE=fopen ( file_ADE , ” w ”) ;
// f p r i n t f ( fp_ADE , ” x ax i s \tWave\n ”) ;
// f p r i n t f ( fp_ADE , ” ( cell )\t ( AU )\n ”) ;
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii++){
f p r i n t f ( fp_ADE , ”%d\ t%l f \n” , i i , Hzin [ nn%2][ i i ] ) ;
}
f c l o s e ( fp_ADE ) ;
}// file ADE writing ends ( i f )
}// f o r ( nn ) *** Time loop ends **********************/
printf (”\ nFinishing Program . . . \ n\n ”) ;
i f (1 ) {// Hz Amplitude output
fp_Hz = fopen (” HzAmp . dat ” , ”w ”) ;
// f p r i n t f ( fp_Hz , ” x−ax i s \ty−ax i s \tHzAmp \n ”) ;
// f p r i n t f ( fp_Hz , ” nm\tnm\tA . U .\ n ”) ;
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii+=2){
f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj+=2){
f p r i n t f ( fp_Hz , ”%l f \ t%l f \ t%l f \n” , xx [ i i ] , yy [ j j ] ,HzAmp[ i i ] [ j j ] ) ;
// column output format
}
}
151
/* f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj++){
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii++)
{
f p r i n t f ( fp_Hz , ”%f \ t ” ,HzAmp[ i i ] [ j j ] ) ;
// Hzshift [ ii ] [ jj ]=Hz [ np0 ] [ ii ] [ jj ] ;
}// matrix output format
f p r i n t f ( fp_Hz , ”\ n ”) ;
}
*/ f c l o s e ( fp_Hz ) ;
}
i f (1 ) {// Hz Time profile ( f o r Fourier DetectorsR4 Analysis ) output , Right Side
fp_Spectrum=fopen (” DetectorsR_Full . dat ” ,” w ”) ;
f p r i n t f ( fp_Spectrum , ” Time ”) ;
f o r ( jj=TFb+1;jj<=metalJ ; jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ” Metal ”) ; // column output format
}
f o r ( jj=metalJ+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ” Air ”) ; // column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
f o r ( int iin=1;iin<=6*(int ) nT ; iin++){// dont change iin++
f p r i n t f ( fp_Spectrum , ”%l f ” , i i n *dt /2/PI ) ;
// 1/(2 PI ) so that the units f o r FFT is rad /attosec
//( omega angular frequency rather than frequency )
f o r ( jj=TFb+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ”%l f ” , DetectorsR4 [ i i n ] [ j j ] ) ;
152
// column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
}
f c l o s e ( fp_Spectrum ) ;
}
i f (1 ) {// Hz Time profile ( f o r Fourier DetectorsR4 Analysis ) output , Right Side
fp_Spectrum=fopen (” DetectorsR_Quarter . dat ” ,” w ”) ;
f p r i n t f ( fp_Spectrum , ” Time ”) ;
f o r ( jj=TFb+1;jj<=metalJ ; jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ” Metal ”) ; // column output format
}
f o r ( jj=metalJ+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ” Air ”) ; // column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
f o r ( int iin=1;iin<=6*(int ) nT ; iin++){
f p r i n t f ( fp_Spectrum , ”%l f ” , i i n *dt /2/PI ) ;
// 1/(2 PI ) so that the units f o r FFT is rad /attosec
//( omega angular frequency rather than frequency )
f o r ( jj=TFb+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ”%l f ” , DetectorsR1 [ i i n ] [ j j ] ) ;
// column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
}
f c l o s e ( fp_Spectrum ) ;
}
153
i f (1 ) {// Hz Time profile ( f o r Fourier DetectorsR1 Analysis ) output , Right Side
fp_Spectrum=fopen (” DetectorsR_Half . dat ” ,” w ”) ;
f p r i n t f ( fp_Spectrum , ” Time ”) ;
f o r ( jj=TFb+1;jj<=metalJ ; jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ” Metal ”) ; // column output format
}
f o r ( jj=metalJ+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ” Air ”) ; // column output format
}
f p r i n t f ( fp_Hz , ”\ n ”) ;
f o r ( int iin=1;iin<=6*(int ) nT ; iin++){
f p r i n t f ( fp_Spectrum , ”%l f ” , i i n *dt /2/PI ) ;
// 1/(2 PI ) so that the units f o r FFT is rad /attosec
//( omega angular frequency rather than frequency )
f o r ( jj=TFb+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ”%l f ” , DetectorsR2 [ i i n ] [ j j ] ) ;
// column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
}
f c l o s e ( fp_Spectrum ) ;
}
i f (1 ) {// Hz Time profile ( f o r Fourier DetectorsR2 Analysis ) output , Right Side
fp_Spectrum=fopen (” DetectorsR_ThreeQuarters . dat ” ,” w ”) ;
f o r ( jj=TFb+1;jj<=metalJ ; jj+=2){
f p r i n t f ( fp_Spectrum , ” Metal\t ”) ; // column output format
}
f o r ( jj=metalJ+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ” Air\t ”) ; // column output format
154
}f p r i n t f ( fp_Hz , ”\ n ”) ;
f o r ( int iin=1;iin<=6*(int ) nT ; iin++){
f p r i n t f ( fp_Spectrum , ”%l f ” , i i n *dt /2/PI ) ;
// 1/(2 PI ) so that the units f o r FFT is rad /attosec
//( omega angular frequency rather than frequency )
f o r ( jj=TFb+1;jj<=TFt−1;jj+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ”%l f ” , DetectorsR3 [ i i n ] [ j j ] ) ;
// column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
}
f c l o s e ( fp_Spectrum ) ;
}
i f (1 ) {// Hz Time profile ( f o r Fourier DetectorsT Analysis ) output , Top Edge
fp_Spectrum=fopen (” DetectorsT . dat ” ,” w ”) ;
f o r ( int iin=1;iin<=6*(int ) nT ; iin++){
f p r i n t f ( fp_Spectrum , ”%l f ” , i i n *dt /2/PI ) ;
// 1/(2 PI ) so that the units f o r FFT is rad /attosec
//( omega angular frequency rather than frequency )
f o r ( ii=centerI ; ii<=SFr ; ii+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ”%l f ” , DetectorsT [ i i n ] [ i i ] ) ;
// column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
}
f c l o s e ( fp_Spectrum ) ;
}
155
i f (1 ) {// Hz Time profile ( f o r Fourier DetectorsB Analysis ) output , Bottom Edge
fp_Spectrum=fopen (” DetectorsB . dat ” ,” w ”) ;
f o r ( int iin=1;iin<=6*(int ) nT ; iin++){
f p r i n t f ( fp_Spectrum , ”%l f ” , i i n *dt /2/PI ) ;
// 1/(2 PI ) so that the units f o r FFT is rad /attosec
//( omega angular frequency rather than frequency )
f o r ( ii=centerI ; ii<=SFr ; ii+=2){
f p r i n t f ( fp_Spectrum , ”\ t ”) ;
f p r i n t f ( fp_Spectrum , ”%l f ” , DetectorsB [ i i n ] [ i i ] ) ;
// column output format
}
f p r i n t f ( fp_Spectrum , ”\ n ”) ;
}
f c l o s e ( fp_Spectrum ) ;
}
i f (1 ) {// Ey Amplitude output
fp_Ey = fopen (” EyAmp . dat ” , ”w ”) ;
// f p r i n t f ( fp_Hz , ” x−ax i s \ty−ax i s \tEyAmp \n ”) ;
// f p r i n t f ( fp_Hz , ” nm\tnm\tA . U .\ n ”) ;
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii+=2){
f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj+=2){
f p r i n t f ( fp_Ey , ”%l f \ t%l f \ t%l f \n” , xx [ i i ] , yy [ j j ] ,EyAmp[ i i ] [ j j ] ) ;
// column output format
}
}
/* f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj++){// Ey Intensity output
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii++)
{
156
f p r i n t f ( fp_Ey , ”%f \ t ” ,EyAmp[ i i ] [ j j ] ) ;
}
f p r i n t f ( fp_Ey , ”\ n ”) ;
}
*/ f c l o s e ( fp_Ey ) ;
}
i f (1 ) {// Ex Amplitude output
fp_Ex = fopen (” ExAmp . dat ” , ”w ”) ;
// f p r i n t f ( fp_Ex , ” x−ax i s \ty−ax i s \tEx \n ”) ;
// f p r i n t f ( fp_Ex , ” nm\tnm\tA . U .\ n ”) ;
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii+=2){
f o r ( jj=Fileoffset ; jj<=SJJ−Fileoffset ; jj+=2){
f p r i n t f ( fp_Ex , ”%l f \ t%l f \ t%l f \n” , xx [ i i ] , yy [ j j ] ,ExAmp[ i i ] [ j j ] ) ;
// column output format
}
}
/* f o r ( jj=0+Fileoffset ; jj<=SJJ−Fileoffset ; jj++){// Ex Intensity output
f o r ( ii=0+Fileoffset ; ii<=SII−Fileoffset ; ii++)
{
f p r i n t f ( fp_Ex , ”%f \ t ” ,ExAmp[ i i ] [ j j ] ) ;
}
f p r i n t f ( fp_Ex , ”\ n ”) ;
}
*/ f c l o s e ( fp_Ex ) ;
}
i f (0 ) {// PML e r r o r test
double error1=0;
double error2=0;
f o r ( ii=SFl+1;ii<=SFr−1;ii++){
error1+=(HzAmp [ ii ] [ centerJ ]−1) *( HzAmp [ ii ] [ centerJ ]−1) ;
157
error2=
error2>abs ( HzAmp [ ii ] [ centerJ ]−1)? error2 : abs ( HzAmp [ ii ] [ centerJ −1]−1) ;
}
cout<<”Total e r r o r=\t”<<s q r t ( error1 )<<endl ;
cout<<”Max e r r o r=\t”<<error2<<endl ;
}
i f (0 ) {// ADE source Amplitude output
fp_Ex = fopen (” HsourceAmp . dat ” , ”w ”) ;
f o r ( ii=SFl+1;ii<=SFr−1;ii++){
f p r i n t f ( fp_Ex , ”%d\ t%l f ” , i i , (HzAmp[ i i ] [ c enterJ ]−1) *100) ;
f p r i n t f ( fp_Ex , ”\ n ”) ;
}
f c l o s e ( fp_Ex ) ;
}
i f (0 ) {// ADE source wave output
fp_Ex = fopen (” Hzin . dat ” , ”w ”) ;
f o r ( jj=CPML_n+Fileoffset ; jj<=SJJ−CPML_n−Fileoffset ; jj++)
{
f p r i n t f ( fp_Ex , ”%d\ t%l f ” , yy [ j j ] , Hzin [ nn%2][ j j ] ) ;
// EHphase [ ii ]=( Hzin [ ( nn+1)%2 ] [ i i ] ) ;//+Hzin [ nn%2][ i i ] ) / 2 . 0 ;
// time average
f p r i n t f ( fp_Ex , ”\ n ”) ;
}
f c l o s e ( fp_Ex ) ;
}
i f (0 ) {// phase i n f o , test
fp_Hz = fopen (” EHphase . dat ” , ”w ”) ;
f p r i n t f ( fp_Hz , ” ii\tEy\ tEyPhase \tHz\tHzPhase \n ”) ;
f p r i n t f ( fp_Hz , ” time\tA . U .\ tDegree \tA . U .\ tDegree \n ”) ;
158
f o r ( ii=Fileoffset ; ii<=SII−Fileoffset ; ii++)
{
}// matrix output format
f p r i n t f ( fp_Hz , ”\ n ”) ;
f c l o s e ( fp_Hz ) ;
}
i f (1 ) {// program i n f o output
fp_info=fopen (” i n f o . txt ” ,” w ”) ;
f p r i n t f ( fp_info , ”\ n−−−−−−−−−−−−−−−−−−−−−−−−−−−Physics \
parameters−−−−−−−−−−−−−−−−−−−−−−−−−−−\n\n ”) ;
f p r i n t f ( fp_info , ” epsilon_0=%f (nF/m) \n” , e p s i l o n 0 ) ;
f p r i n t f ( fp_info , ” mu_0=%f (nH/m) \n” ,mu 0) ;
f p r i n t f ( fp_info , ” c_0=%f (nm/ at to s e c ) \n” , c 0 ) ;
f p r i n t f ( fp_info , ” lambda_0=%f (nm) \n” , lambda 0 ) ;
f p r i n t f ( fp_info , ” omega_0=%f ( rad / a t to s e c ) \n” , omega 0 ) ;
f p r i n t f ( fp_info , ” Frequency f_0=%e (GGHz) \n” , omega 0 /2 ./PI ) ;
f p r i n t f ( fp_info , ” k_0=%f (/nm) \n” , k 0 ) ;
f p r i n t f ( fp_info , ” epsilon_inf=%f \n\n” , e p s i l o n i n f ) ;
f p r i n t f ( fp_info , ” ( Debye )\ tdepsilon_p=%f \n” , d ep s i l on p ) ;
f p r i n t f ( fp_info , ” ( Debye )\ttau_p=%f \n” , tau p ) ;
f p r i n t f ( fp_info , ” ( Lorentz )\tdelta_p=%f \n” , d e l t a p ) ;
f p r i n t f ( fp_info , ” ( Drude )\ tomega_p=%f \n” , omega p ) ;
f p r i n t f ( fp_info , ” ( Drude )\ tgamma_p=%f \n” ,gamma p) ;
f p r i n t f ( fp_info , ” ( Kerr )\talpha=%f \n” , alpha ) ;
f p r i n t f ( fp_info , ” ( Kerr )\tchi (3)=%e (m2/V2) \n” , chi3M) ;
f p r i n t f ( fp_info , ” ( Raman )\ttau1=%f \n” , tau1 ) ;
f p r i n t f ( fp_info , ” ( Raman )\ttau2=%f \n” , tau2 ) ;
f p r i n t f ( fp_info , ” ( Raman )\ tomega_raman=%f \n” , omega raman ) ;
f p r i n t f ( fp_info , ” ( Raman )\ tdelta_raman=%f \n\n” , del ta raman ) ;
f p r i n t f ( fp_info , ” For metal , Re ( epsilon )=%f \tIm ( ep s i l o n )=%f \
159
\n ” , epsilon_R , epsilon_I ) ;
f p r i n t f ( fp_info , ” Hz input amplitude :\ t%E (A/m) \n” , smoothupHc ) ;
f p r i n t f ( fp_info , ” SPP wave l ength in medium lambda=%f (nm) \n” , lambda) ;
f p r i n t f ( fp_info , ” SPP wave vector Re ( kx )=%f (1/nm) \tIm (kx )=%e (1/nm) \n\
\n ” , kx_R , kx_I ) ;
f p r i n t f ( fp_info , ” Evanescece wave in air \t Re ( ky0 )=%f (1/nm) \ t \tIm ( ky0 )=%e\
(1/ nm ) \n Or penetration depth\t 1/Re ( ky0 )=%f (nm) \ t \ t1 /Im\
( ky0 )=%f (nm) \n” , ky0 R , ky0 I , 1 . / ky0 R , 1 . / ky0 I ) ;
f p r i n t f ( fp_info , ”\ nEvanescece wave in metal\t \
Re ( ky2 )=%f (1/nm) \ t \tIm ( ky2 )=%e (1/nm) \n \
Or penetration depth \t 1/Re ( ky2 )=%f (nm) \ t \ t1 /Im( ky2 )=\
%f (nm) \n” , ky2 R , ky2 I , 1 . / ky2 R , 1 . / ky2 I ) ;
f p r i n t f ( fp_info , ” Ey input amplitude (V/m ) :\ t ( in air ) %E\
\t ( in metal ) %E\n” ,\
smoothupEc *abs_kx /omega_0 /epsilon_0 /epsilon_inf ,
smoothupEc *abs_kx /omega_0 /epsilon_0 ) ;
f p r i n t f ( fp_info , ” Incident wave period :%d time s t ep s (=%f f s ) \n\n” ,
int ( nT ) , nT*dt /1000) ;
f p r i n t f ( fp_info , ”\ n\n\n−−−−−−−−−−−−−−−−−−−−−−−−−−−FDTD \
parameters−−−−−−−−−−−−−−−−−−−−−−−−−−−\n\n ”) ;
f p r i n t f ( fp_info ,
”Total time simulated : steps=%d (=%f f s ) \n” , Ntotal , dt *Ntotal /1000 . ) ;
f p r i n t f ( fp_info , ” dt=%f attosecond\n” , dt ) ;
f p r i n t f ( fp_info , ” File written every %d time s t ep s (=%f f s ) \
\twith ( ABC frame ) offset=%d\n\n\n” ,
filestep , dt*filestep /1000 . , Fileoffset ) ;
f p r i n t f ( fp_info , ” Total cell numbers along x : %d\ t \ ty : %d\n” , SII , SJJ ) ;
f p r i n t f ( fp_info , ” Grids ( coarse ) sampling density : ( x )%d\ t ( y )%d\n” ,
160
Nlambdax , Nlambday ) ;
f p r i n t f ( fp_info , ” dxmax=%f nm\ tdxmin=%f nm\tx−(d i s t ance ) t o t a l=%f nm\n” ,
dxmax , dxmin , xx [ SII ] ) ;
f p r i n t f ( fp_info , ” dymax=%f nm\ tdymin=%f nm\ty−(d i s t ance ) t o t a l=%f nm\n” ,
dymax , dymin , yy [ SJJ ] ) ;
f p r i n t f ( fp_info , ” Air−Metal interface at y ( J )=%d(+) ˜ %f (nm) \n” ,
metalJ , yy [ metalJ ] ) ;
f p r i n t f ( fp_info , ” Total field frame :\ t ( Left , Right , Bottom , Top )=\
(%d,%d,%d,%d) ( c e l l ) ˜(% f ,%f ,%f ,% f ) (nm) \n” ,
TFl , TFr , TFb , TFt , xx [ TFl ] , xx [ TFr ] , yy [ TFb ] , yy [ TFt ] ) ;
f p r i n t f ( fp_info , ” Scattered field frame :\ t \
( Left , Right , Bottom , Top )=\
(%d,%d,%d,%d) ( c e l l ) ˜(% f ,%f ,%f ,% f ) (nm) \n” ,
SFl , SFr , SFb , SFt , xx [ SFl ] , xx [ SFr ] , yy [ SFb ] , yy [ SFt ] ) ;
f p r i n t f ( fp_info , ” CPML layer width CPML_n=%d\n” ,CPML n) ;
f p r i n t f ( fp_info , ”\ nTotal nonuniform s i z e numbers :\
x %d\ t y %d\n” ,0 , varygrade n ) ;
f p r i n t f ( fp_info , ” s i z e ratio of max to min cell maxVSmin=%f \n” ,maxVSmin) ;
f p r i n t f ( fp_info , ” s i z e ratio of adjacent cells : %f or %f \n” ,
gradeRatey , 1 . 0 / gradeRatey ) ;
f p r i n t f ( fp_info , ” nonuniform area ( y ) : %d\\%d %d/%d\
( cell ) \t˜\t%f\\% f %f/%f (nm) \n” ,
levely0 , levely1 , levely2 , levely3 , yy [ levely0 ] ,
yy [ levely1 ] , yy [ levely2 ] , yy [ levely3 ] ) ;
f p r i n t f ( fp_info , ” Field Ramp Up %d ˜wavelength ( s ) \n\n” ,
( int ) s q r t ( ( double ) rampup ) ) ;
f p r i n t f ( fp_info , ”\ n\ndebye model FDTD parameters : a=%f \ tb=%f \
\n ” , a_debye , b_debye ) ;
f p r i n t f ( fp_info , ” lorentz model FDTD parameters : a=%f \ tb=%f \ tc=%f \
\n ” , a_lorentz , b_lorentz , c_lorentz ) ;
f p r i n t f ( fp_info , ” drude model FDTD parameters : a=%f \ tb=%f \ tc=%f \
161
\n ” , a_drude , b_drude , c_drude ) ;
f p r i n t f ( fp_info , ” Raman model FDTD parameters : a=%f \ tb=%f \ tc=%f \n\
\n ” , a_raman , b_raman , c_raman ) ;
f p r i n t f ( fp_info , ” CPML factors : sigma_max=%f \ tkkx max=%f \taax max=\
%f \n\n” , sigma max , kkx max , aax max ) ;
f p r i n t f ( fp_info , ” Nonlinear iteration loop number : nonlinearM=\
%d\n\n” , nonl inearM) ;
f p r i n t f ( fp_info , ”\ n\n\n−−−−−−−−−−−−−−−−−−−−−−−−−−−CPML \
parameters−−−−−−−−−−−−−−−−−−−−−−−−−−−\n\n ”) ;
f p r i n t f ( fp_info , ” sigma_max=%f \n” , sigma max ) ;
f p r i n t f ( fp_info , ” kkx_max=%f \n” , kkx max ) ;
f p r i n t f ( fp_info , ” aax_max=%f \n” , aax max ) ;
f p r i n t f ( fp_info , ”\ n\n\n−−−−−−−−−−−−−−−−−−−−−−−−−−−Time/Date \
i n f o−−−−−−−−−−−−−−−−−−−−−−−−−−−\n\n ”) ;
time_t rawtime ;
struct tm * timeinfo ;
time ( &rawtime ) ;
timeinfo = localtime ( &rawtime ) ;
f p r i n t f ( fp_info , ”\ nThis file was written at ( date /time ) %s ” ,
asctime ( timeinfo ) ) ;
f c l o s e ( fp_info ) ;
}
printf (”\ t\t\tFile written f o r %s \n” ,” Average I n t e n s i t y ”) ;
// step indicator on screen
/*−−−−Endline of . dat file writing ***************************/
162
/*−−−−−−−−−−−−−−−−−−−−−−−−clean up−−−−−−−−−−−−−−−−−−−−−−−*/
free3D (Hz , 3 , SII+1) ;//1
free3D ( Hzx , 3 , SII+1) ;//2
free3D ( Hzy , 3 , SII+1) ;//3
free3D (Dy , 3 , SII+1) ; / / 2 . 2 Displacement y
free3D (Dx , 3 , SII+1) ; / / 3 . 2 Displacement x
free3D ( ET2 , 3 , SII+1) ; / / 3 . 4 | E |ˆ2
free3D (Ey , 3 , SII+1) ;//4
free3D (Ex , 3 , SII+1) ;//5
free3D ( Px_Debye , 3 , SII+1) ;
free3D ( Px_Lorentz , 3 , SII+1) ;
free3D ( Px_Raman , 3 , SII+1) ;
free3D ( Px_Kerr , 3 , SII+1) ;
free3D ( Px_Drude , 3 , SII+1) ;
free3D ( Py_Drude , 3 , SII+1) ;
free3D ( Py_Debye , 3 , SII+1) ;
free3D ( Py_Lorentz , 3 , SII+1) ;
free3D ( Py_Kerr , 3 , SII+1) ;
free3D ( Py_Raman , 3 , SII+1) ;
free3D (Sp , 3 , SII+1) ;
// cout<<”11”<<endl ;
free2D ( epsilon , SII+1) ;//6
free1D ( sigmax ) ; //8
free1D ( sigmay ) ; //9
free2D (mu , SII+1) ;//7
// free2D ( sigmastarx , SII+1) ;//10
free3D ( psi_dyHz , 3 , SII+1) ;//11
free3D ( psi_dxHz , 3 , SII+1) ;//
free3D ( psi_dxEy , 3 , SII+1) ;//
free3D ( psi_dyEx , 3 , SII+1) ;//
free1D ( aax ) ; free1D ( bbx ) ; free1D ( ccx ) ;
free1D ( aay ) ; free1D ( bby ) ; free1D ( ccy ) ;
163
// free1D ( a_Hz ) ; free1D ( b_Hz ) ; free1D ( c_Hz ) ;
free1D ( kkx ) ; free1D ( kky ) ; // free1D ( k_Hz ) ;
free2D ( Hzin , 2 ) ;//18
free2D ( Eyin , 2 ) ;//19
free2D ( HzAmp , SII+1) ;//20
free2D ( ExAmp , SII+1) ;//21
free2D ( EyAmp , SII+1) ;//22
free2D ( DetectorsR4 , 6* ( int ) nT+2) ;
free2D ( DetectorsR1 , 6* ( int ) nT+2) ;
free2D ( DetectorsR2 , 6* ( int ) nT+2) ;
free2D ( DetectorsR3 , 6* ( int ) nT+2) ;
free2D ( DetectorsT , 6* ( int ) nT+2) ;
free2D ( DetectorsB , 6* ( int ) nT+2) ;
free1D ( HsourceAmp ) ;//23
free1D ( EHphase ) ;//24
//
free1D ( dx ) ;//38
free1D ( dy ) ;//39
free1D ( hxi ) ;//40
free1D ( hyj ) ;//41
free1D ( xx ) ;
free1D ( yy ) ;
printf (”\ n ”) ;
printf (” File written successfully . Now closing it . . . \ n ”) ;
printf (”\ n ”) ;
// system (” PAUSE ”) ;
r e turn EXIT_SUCCESS ; / / ( main )
}// end main ( )
/*** Endline of main program body **************************************/
164
double Ex_pml_func ( double Ex_previous , double Hzup , double Hzdown , double sigma ,
double epsilon , double d_x ) {// general E ( n+1)calc in CPML_n layer
double CE=(epsilon−(sigma*dt *0.5/ epsilon_0 ) ) /( epsilon+
( sigma *dt *0.5/ epsilon_0 ) ) ; // factor f o r E
double CH=(dt/ epsilon_0 ) /( epsilon+sigma*dt *0.5/ epsilon_0 ) /d_x ; // factor f o r H
double Ex=CE*Ex_previous+CH *( Hzup−Hzdown ) ;
//H_2 has larger space index than H_1
r e turn Ex ;
}// Ex_pml_func
double Ey_pml_func ( double Ey_previous , double Hzright , double Hzleft ,
double sigma , double epsilon , double d_y ) {// general E (n+1)calc in CPML_n layer
double CE=(epsilon−(sigma*dt *0.5/ epsilon_0 ) ) /( epsilon+
( sigma *dt *0.5/ epsilon_0 ) ) ;
double CH=−(dt/epsilon_0 ) /( epsilon+sigma *dt *0.5/ epsilon_0 ) /d_y ;
double Ey=CE*Ey_previous+CH *( Hzright−Hzleft ) ;
r e turn Ey ;
}// Ey_pml_func
double Hzx_func ( double H_previous , double Eyright , double Eyleft , double sigmastar ,
double mu , double d_x ) {// gneral H ( n+1) calc in CPML_n layer
double CH=(mu−sigmastar *dt *0.5/ mu_0 ) /(1+ sigmastar *dt *0.5/ mu_0 ) ;
double CE=−(dt/mu_0 ) /( mu+sigmastar *dt *0.5/ mu_0 ) /d_x ;
double Hz=CH*H_previous+CE *( Eyright−Eyleft ) ;
r e turn Hz ;
}// Hzx_func
double Hzy_func ( double H_previous , double Exup , double Exdown , double sigmastar ,
double mu , double d_y ) {// gneral H ( n+1) calc in CPML_n layer
double CH=(mu−sigmastar *dt *0.5/ mu_0 ) /( mu+sigmastar *dt *0.5/ mu_0 ) ;
double CE=(dt/mu_0 ) /( mu+sigmastar *dt *0.5/ mu_0 ) /d_y ;
165
double Hz=CH*H_previous+CE *( Exup−Exdown ) ;
r e turn Hz ;
}// Hzy_func
void Hz_contour_func ( double *** Hz , double *** Ex , double *** Ey , int time , int ii ,
int jj , double area , double Fn , double Fs , double Fw , double Fe ){
//North , South , West , East ; Integral Form of Maxwell Equations , contour calc
int pretime=(time+2)%3 ;
double En=Ex [ pretime ] [ ii ] [ jj ] ;
double Es=Ex [ pretime ] [ ii ] [ jj−1] ;
double Ew=Ey [ pretime ] [ ii−1] [ jj ] ;
double Ee=Ey [ pretime ] [ ii ] [ jj ] ;
Hz [ time ] [ ii ] [ jj ]=Hz [ pretime ] [ ii ] [ jj ]+dt/mu_0/area * ( ( En*Fn−Es*Fs ) *dxmax+
( Ew*Fw−Ee*Fe ) *dymax ) ;
}// Hz_contour_func
int SFflag ( int ii , int jj ){
int f l a g=(ii<=SFr )&&(ii>=SFl )&&(jj<=SFt )&&(jj>=SFb ) ;
// Boolean f l a g , 1 f o r in Scattered fields area , 0 f o r NOT .
r e turn f l a g ;
}
int TFflag ( int ii , int jj ){
int f l a g=(ii<=TFr )&&(ii>=TFl )&&(jj<=TFt )&&(jj>=TFb ) ;
// Boolean f l a g , 1 f o r in Total fields area , 0 f o r NOT .
r e turn f l a g ;
}
int CPMLflag ( int ii , int jj ) {
int f l a g=(ii<SFl )&&(ii>SFr )&&(jj<SFb )&&(jj>SFt ) ;
// Boolean f l a g , 1 f o r in CPML area , 0 f o r NOT .
r e turn f l a g ;
166
}int Metalflag ( int ii , int jj ) {// determine i f (ii , jj ) is in metal
int f l a g=(jj<=metalJ ) ; / / ( ii<=SFr )&&(ii>=SFl )&&(jj<=metalJ )&&(jj>=SFb ) ;
r e turn f l a g ;
}
int Airflag ( int ii , int jj ) {// determine i f (ii , jj ) is in air
int f l a g=(jj>metalJ ) ; / / ( ii<=SFr )&&(ii>=SFl )&&(jj>metalJ )&&(jj<=SFt ) ;
r e turn f l a g ;
}
FDTDarray.h
#ifndef FDTDarray_h
#define FDTDarray_h
double * array1D ( int n ) ;
double ** array2D ( int narrows , int ncolumns ) ;
double *** array3D ( int narrows , int ncolumns , int nz ) ;
void free1D ( double * p1 ) ;
void free2D ( double ** p2 , int nx ) ;
void free3D ( double *** p3 , int nx , int ny ) ;
double array_max ( double * p , int from , int to ) ;
double array_min ( double * p , int from , int to ) ;
double array_amp ( double * p , int from , int to ) ;
167
double Dx_func ( double Dx_previous , double Hzxdown , double Hzxup , double hyj , double dt ) ;
double Dy_func ( double Dy_previous , double Hzyleft , double Hzyright ,
double hxi , double dt ) ;
double sqrtReal ( double a , double b , double c , double d ) ;
double sqrtImag ( double a , double b , double c , double d ) ;
double sqrtReal ( double a , double b ) ;
double sqrtImag ( double a , double b ) ;
double square ( double a ) ;
double cube ( double a ) ;
#endif
168
FDTDarray.cpp
/* Malloc_multi_array */
/* malloc to defined a 2 or 3 dimensional array */
#include <stdio . h>
#include <iostream>
#include <stdlib . h>
#include <cmath>
#include ”FDTDarray . h”
using namespace std ;
/*−−−−−−−−−−−−−−end #include−−−−−−−−−−*/
/*−−−−−−−−−−−−−−start array1D f unc t i on definition−−−−−−*/
double * array1D ( int n )
{
int i=0;
double * p1 ;
p1=(double *) malloc (n*sizeof ( double ) ) ;
i f ( p1==NULL ){
f p r i n t f ( stderr , ” out of memory \n ”) ;
r e turn 0 ;
}// i f
f o r ( i=0;i<=n−1;i++){
p1 [ i ]=0;
}// f o r ( i )
r e turn ( p1 ) ;
}// array1D ( )
169
/*−−−−−−−−−−−−−−end array1D f unc t i on definition−−−−−−*/
/*−−−−−−−−−−−−−−−−−−−Start array2D f unc t i on definition−*/
double ** array2D ( int nrows , int ncolumns )
{
int i=0;
int j=0;
double **p2 ;
p2 = ( double **) malloc ( nrows * sizeof ( double *) ) ;
i f ( p2 == NULL ) {
f p r i n t f ( stderr , ”out of memory \n ”) ;
r e turn 0 ;
}//( i f )
f o r ( i = 0 ; i <= nrows−1; i++){
p2 [ i ] = ( double *) malloc ( ncolumns * sizeof ( double ) ) ;
i f ( p2 [ i ] == NULL ){
f p r i n t f ( stderr , ”out of memory \n ”) ;
r e turn 0 ;
}//( i f )
}// ( f o r ( i ) )
f o r ( i=0;i<=nrows−1;i++){
f o r (j=0;j<=ncolumns−1;j++){
p2 [ i ] [ j ]=0;
}//( f o r ( i ) )
}//( f o r ( j ) )
r e turn ( p2 ) ;
}// array2D
/*−−−−−−−−−−−−−−−−−−−end array2D definitino−−−−−−−−−−−−−−−−*/
170
/*−−−−−−−−−−−−−−−−−−−Start array3D f unc t i on definition−*/
double *** array3D ( int nrows , int ncolumns , int nz )
{
int i=0;
int j=0;
int k=0;
double *** p3 ;
p3 = ( double ***) malloc ( nrows * sizeof ( double **) ) ;
i f ( p3 == NULL ) {
f p r i n t f ( stderr , ”out of memory \n ”) ;
r e turn 0 ;
}//( i f )
f o r ( i = 0 ; i <= nrows−1; i++){
p3 [ i ] = ( double **) malloc ( ncolumns * sizeof ( double *) ) ;
i f ( p3 [ i ] == NULL ){
f p r i n t f ( stderr , ”out of memory \n ”) ;
r e turn 0 ;
}//( i f )
}
f o r ( i = 0 ; i <= nrows−1; i++)
f o r ( int j=0;j<=ncolumns−1;j++){
p3 [ i ] [ j ] = ( double *) malloc ( nz * sizeof ( double ) ) ;
i f ( p3 [ i ] [ j ] == NULL ){
f p r i n t f ( stderr , ”out of memory \n ”) ;
r e turn 0 ;
}//( i f )
}// ( f o r ( i ) )
171
f o r ( j=0;j<=nrows−1;j++)
f o r (i=0;i<=ncolumns−1;i++)
f o r (k=0;k<=nz−1;k++){
p3 [ j ] [ i ] [ k ]=0;
}
r e turn ( p3 ) ;
}// array3D
/*−−−−−−−−−−−−−−−−−−−end array3D definitino−−−−−−−−−−−−−−−−*/
/*−−−−−−−−−−−−−−−−−−free alloc memory−−−−−−−−−−−−−−−−−−*/
void free1D ( double * p1 ){
free ( p1 ) ;
// printf (”\ nMemory is freed . . . \ n ”) ;
}
void free2D ( double ** p2 , int nx ){
f o r ( int i=0;i<=nx−1;i++){
free ( p2 [ i ] ) ;
}
free ( p2 ) ;
// printf (”\ nMemory is freed . . . \ n ”) ;
}
void free3D ( double *** p3 , int nx , int ny ){
f o r ( int i=0;i<=nx−1;i++)
f o r ( int j=0;j<=ny−1;j++)
free ( p3 [ i ] [ j ] ) ;
f o r ( int k=0;k<=nx−1;k++)
free ( p3 [ k ] ) ;
172
// printf (”\ nMemory is freed . . . \ n ”) ;
free ( p3 ) ;
}
/*−−−−−−−−−−−−−−−−−−−array member operation−−−−−−−−−−−−−−*/
double array_max ( double * p , int from , int to ){
double maxx=0;
f o r ( int i=from ; i<=to ; i++){
maxx=maxx>p [ i ] ? maxx : p [ i ] ;
}
r e turn maxx ;
}
double array_min ( double * p , int from , int to ){
double minn=0;
f o r ( int i=from ; i<=to ; i++){
minn=minn<p [ i ] ? minn : p [ i ] ;
}
r e turn minn ;
}
double array_amp ( double * p , int from , int to ){
double amp=0.0;
f o r ( int i=from ; i<=to ; i++){
amp=amp>fabs ( p [ i ] ) ?amp : fabs ( p [ i ] ) ;
}
r e turn amp ;
}
/*−−−−−−−−−−−FDTD Displacement vector functions−−−−−−−−−−−−−−−*/
double Dx_func ( double Dx_previous , double Hzxdown , double Hzxup , double hyj ,
173
double dt ) {
double Dx=0;
Dx=Dx_previous+dt/hyj *( Hzxup−Hzxdown ) ;
r e turn Dx ;
}
double Dy_func ( double Dy_previous , double Hzyleft , double Hzyright , double hxi ,
double dt ) {
double Dy=0;
Dy=Dy_previous+dt/hxi *( Hzyleft−Hzyright ) ;
r e turn Dy ;
}
/*−−−−−−−−−−−−−−caculate Real and Imag parts−−−−−−−−−−−−−−−−*/
/* assume (x+yi ) ˆ2=(a+bi ) /(c+di ) , and x , y unknow whi le a , b , c , d
given−positive root only , which means −PI<arg (z )<=PI−−*/
double sqrtReal ( double a , double b , double c , double d ){
double e=(a*c+b*d ) /( c*c+d*d ) ; // r e a l part after squared
double f=(b*c−a*d ) /( c*c+d*d ) ; // imag part after squared
r e turn sq r t ( 0 . 5 * ( e+sqr t ( e*e+f*f ) ) ) ;
}
double sqrtImag ( double a , double b , double c , double d ){
double e=(a*c+b*d ) /( c*c+d*d ) ; // r e a l part after squared
double f=(b*c−a*d ) /( c*c+d*d ) ; // imag part after squared
r e turn sq r t (0.5*(−e+sqr t (e*e+f*f ) ) ) ;
}
double sqrtReal ( double a , double b ){
r e turn sq r t ( ( a+sqr t ( a*a+b*b ) ) *0 . 5 ) ;
}
double sqrtImag ( double a , double b ){
r e turn sq r t ((−a+sqr t ( a*a+b*b ) ) *0 . 5 ) ;
}
174
175
FDTDcubic.h
#ifndef FDTDcubic_h
#define FDTDcubic_h
#include <complex>
using namespace std ;
double square ( double a ) ;
double cube ( double a ) ;
complex<double> cbrt ( complex<double>) ; // reload cubic root f unc t i on
complex<double> cartesian ( double x , double y ) ; // complex number assign
complex<double> cartesian ( complex<double> z ) ;
class LinearFormula{// generic Linear Equation
double _k , _b , _root ; // kx+b=0;
public :
LinearFormula ( ) ; // default constructor
LinearFormula ( double , double ) ; // constructor
int setFormula ( double , double ) ; // change _k , _b i f necessary
int realRootNumber ( ) ;
// f i nd how many r e a l r oo t s exit , reload f o r 1st , 2 nd and 3rd order
double linearRoot ( ) ; // r eturn solution
} ;
class QuadraticFormula {// generica Quadratic Equation , axˆ2+bx+c==0 with a=1;
double _b , _c ; // xˆ2+bx+c==0,generic form f o r quadratic functions
double _delta ; // discriminant , = bˆ2−4ac f o r a=1
complex<double> _root1 , _root2 ;
public :
176
QuadraticFormula ( ) ; // default constructor
QuadraticFormula ( double , double ) ; // constructor
int setFormula ( double , double ) ; // change , _b , _c , _delta
int realRootNumber ( ) ;
complex<double> quadraticRoot1 ( ) ;
complex<double> quadraticRoot2 ( ) ;
} ;
class FullQuadraticFormula {// general Quadratic Equation , axˆ2+bx+c=0
//with the option that a==0;
double a_ , b_ , c_ ; // xˆ2+bx+c==0,generic form f o r quadratic functions
double delta_ ; // discriminant , = bˆ2−4ac f o r a=1
int isQuadratic_ ; // f l a g to show i f a !=0( isQuadratic=1,true )
complex<double> root1_ , root2_ ;
public :
FullQuadraticFormula ( ) ; // default constructor
FullQuadraticFormula ( double , double , double ) ; // constructor
int setFormula ( double , double , double ) ;
//change , a_ , b_ , c_ , delta_
int realRootNumber ( ) ;
complex<double> quadraticRoot1 ( ) ;
complex<double> quadraticRoot2 ( ) ;
} ;
class CubicFormula{
protected :
double _p , _q ; // xˆ3+ px+ q=0, canonic form
double _discriminant ; // discriminant , = (p /3) ˆ3+(q /2) ˆ2 ,
//d<0 : 3 r e a l r oo t s ;
//d==0: 3 r e a l r oo t s and at least 2 are equal
//d>0 : 1 r e a l root and 2 complex r oo t s
complex<double> _sqrtD ;
177
// square root of the discriminant , which is the source f o r complex r oo t s
complex<double> _uu , _vv ; // intemedium variables , part of _u and _v
complex<double> _u , _v ; // intemedium variables , x=u+v
complex<double> CROOT2 ; // prefactor , −1/2+i* s q r t (3) /2
complex<double> CROOT3 ;
public :
CubicFormula ( ) ; // default constructor
CubicFormula ( double , double ) ;
// constructor to initialize canonic cubic formula
int setFormula ( double , double ) ;
int realRootNumber ( ) ;
complex<double> cubicRoot1 ( ) ;
complex<double> cubicRoot2 ( ) ;
complex<double> cubicRoot3 ( ) ;
} ;
/*Full cubicFormula*/
class FullCubicFormula {// a0*xˆ3+axˆ2+bx+c=0; general form
protected :
double a0_ , a_ ; //
double p_ , q_ ; // xˆ3+ px+ q=0, canonic form
double discriminant_ ; // discriminant , = (p /3) ˆ3+(q /2) ˆ2 ,
double isCubic_ ; // f l a g a0 !=0 ( true , 1)
//d<0 : 3 r e a l r oo t s ;
//d==0: 3 r e a l r oo t s and at least 2 are equal
//d>0 : 1 r e a l root and 2 complex r oo t s
complex<double> sqrtD_ ;
// square root of the discriminant , which is the source f o r complex r oo t s
complex<double> uu_ , vv_ ; // intemedium variables , part of _u and _v
complex<double> u_ , v_ ; // intemedium variables , x=u+v
complex<double> CROOT2 ;
complex<double> CROOT3 ;
178
complex<double> root1_ ;
complex<double> root2_ ;
complex<double> root3_ ;
public :
FullCubicFormula ( ) ; // default constructor
FullCubicFormula ( double , double , double , double ) ;
// constructor to initialize canonic cubic formula
int setFormula ( double , double , double , double ) ;
int realRootNumber ( ) ;
complex<double> cubicRoot1 ( ) ;
complex<double> cubicRoot2 ( ) ;
complex<double> cubicRoot3 ( ) ;
} ;
#endif
179
FDTDcubic.cpp
#include <cmath>
#include <iostream>
#include <complex>
#include ”FDTDcubic . h”
using namespace std ;
double square ( double a ){
r e turn a*a ;
}
double cube ( double a ) {
r e turn a*a*a ;
}
complex<double> cbrt ( complex<double> z ) {
r e turn pow (z , 1 . 0 / 3 ) ;
}
/* Cartesian Cast Function , overloaded */
complex<double> cartesian ( double x , double y ) {
r e turn complex<double> (x , y ) ;
}
complex<double> cartesian ( complex<double> z ) {
r e turn complex<double> ( r e a l ( z ) , imag ( z ) ) ;
}
//kx+b==0;
/* linearFormular Class , to get the slope ( root ) ************************/
LinearFormula : : LinearFormula ( ) {// default constructor , pascal ' s triangle
180
_k=1;
_b=1;
// default solution = −1
}
LinearFormula : : LinearFormula ( double k , double b ) {
setFormula (k , b ) ;
}
int LinearFormula : : realRootNumber ( ) {
i f ( _k !=0){
_root=−_b/_k ;
r e turn 1;//1 root
}
e l s e {
r e turn 0 ; // no root
}
}
double LinearFormula : : linearRoot ( ) {
i f ( _k !=0){
r e turn (−_b/_k ) ;
}
e l s e {
r e turn (−1) ;
}
}
int LinearFormula : : setFormula ( double k , double b ){
_k=k ;
_b=b ;
i f (_k>0)
r eturn 1 ;
e l s e i f ( _k<0)
r eturn −1;
e l s e {
181
cout<<”Warning : Slope is zero , no f i n i t e root !\ n”<<endl ;
r e turn 0 ;
}
}
//xˆ2+bx+c==0;
/* QuadraticFormula Class , 2 roots , double r oo t s or complex r oo t s *************/
QuadraticFormula : : QuadraticFormula ( ) {// default constructor
// using Pascal ' s triangle as coefficients , which gives two identical r oo t s −1;
setFormula ( 2 , 1 ) ;
}
QuadraticFormula : : QuadraticFormula ( double b , double c ) {
setFormula (b , c ) ;
}
int QuadraticFormula : : realRootNumber ( ) {// number of r e a l r oo t s
i f ( _delta>0)
r eturn 2 ;
e l s e i f ( _delta<0)
r eturn 0 ;
e l s e
r e turn 1 ;
}
int QuadraticFormula : : setFormula ( double b , double c ) {//xˆ2+bx+c==0
_b=b ;
_c=c ;
_delta=_b*_b−4.0* _c ;
i f ( _delta>0)
r eturn 1 ;
e l s e i f ( _delta<0)
r eturn −1;
e l s e
r e turn 0 ;
}
182
complex<double> QuadraticFormula : : quadraticRoot1 ( ) {
i f ( _delta>=0)
_root1=(−_b+sqr t ( _delta ) ) / 2 . 0 ;
e l s e
_root1=cartesian (−_b /2 . 0 , s q r t (−_delta ) /2 . 0 ) ;
r e turn ( _root1 ) ;
}
complex<double> QuadraticFormula : : quadraticRoot2 ( ) {
i f ( _delta>=0)
_root2=(−_b−s q r t ( _delta ) ) / 2 . 0 ;
e l s e
_root2=cartesian (−_b /2.0 ,− s q r t (−_delta ) /2 . 0 ) ;
r e turn ( _root2 ) ;
}
//axˆ2+bx+c==0;
/* FullQuadraticFormula Class , 2 roots , double r oo t s or complex r oo t s *************/
/* taking account f o r the case a=0*/
FullQuadraticFormula : : FullQuadraticFormula ( ) {// default constructor
// using Pascal ' s triangle as coefficients , which gives two identical r oo t s −1;
setFormula ( 1 , 2 , 1 ) ;
}
FullQuadraticFormula : : FullQuadraticFormula ( double a , double b , double c ){
setFormula (a , b , c ) ;
}
int FullQuadraticFormula : : realRootNumber ( ) {// number of r e a l r oo t s
i f ( delta_>0)
r eturn 2 ;
e l s e i f ( delta_<0)
r eturn 0 ;
e l s e
r e turn 1 ;
}
183
int FullQuadraticFormula : : setFormula ( double a , double b , double c ) {// axˆ2+bx+c==0
isQuadratic_=1;
a_=a ;
i f (a==0){
// cout<<”Warning : Quadratic term vanishes !”<<endl ;
LinearFormula linearTemp (b , c ) ;
root1_=linearTemp . linearRoot ( ) ;
root2_=root1_ ;
isQuadratic_=0;
r eturn −1;
}
e l s e {
b_=b/a ;
c_=c/a ;
delta_=b_*b_−4.0* c_ ;
i f ( delta_>=0)
r eturn 1 ;
e l s e
r e turn 0 ;
}
}
complex<double> FullQuadraticFormula : : quadraticRoot1 ( ) {
i f ( isQuadratic_ !=0){
i f ( delta_>=0)
root1_=(−b_+sqr t ( delta_ ) ) / 2 . 0 ;
e l s e
root1_=cartesian (−b_ /2 . 0 , s q r t (−delta_ ) /2 . 0 ) ;
}
r e turn ( root1_ ) ;
}
complex<double> FullQuadraticFormula : : quadraticRoot2 ( ) {
184
i f ( isQuadratic_ !=0){
i f ( delta_>=0)
root2_=(−b_−s q r t ( delta_ ) ) / 2 . 0 ;
e l s e
root2_=cartesian (−b_ /2.0 ,− s q r t (−delta_ ) /2 . 0 ) ;
}
r e turn ( root2_ ) ;
}
/* CubicFormula Class , 3 r oo t s in total ***********************************/
CubicFormula : : CubicFormula ( double p , double q ){
CROOT2=cartesian (−0.5 ,0 .8660254037844386) ;//−1/2+i* s q r t (3) /2
CROOT3=cartesian (−0.5 ,−0.8660254037844386) ;//−1/2−i* s q r t (3) /2
setFormula (p , q ) ;
}
int CubicFormula : : setFormula ( double p , double q ) {
_p=p ;
_q=q ;
int returnValue =0;
_discriminant=cube ( _p /3 . 0 )+square ( _q /2 . 0 ) ;
// http : // en . wikipedia . org/wiki/ Cubic_equation
i f ( _discriminant==0){
_sqrtD=0;
_uu=cbrt(−_q /2 . 0 ) ;
_vv=cbrt(−_q /2 . 0 ) ;
returnValue =0;
}
i f ( _discriminant >0){
_sqrtD=sqr t ( _discriminant ) ;
_uu=cbrt ( r e a l (−_q/2.0+ _sqrtD ) ) ;
_vv=cbrt ( r e a l (−_q/2.0− _sqrtD ) ) ;
185
returnValue =1;
}
e l s e {
_sqrtD=cartesian (0 , s q r t (−_discriminant ) ) ;
_uu=cbrt(−_q/2.0+ _sqrtD ) ;
_vv=cbrt(−_q/2.0− _sqrtD ) ;
returnValue=−1;
}
r e turn returnValue ;
}
int CubicFormula : : realRootNumber ( ) {
i f ( _discriminant <0)
r eturn 3 ;
e l s e
r e turn 1 ;
}
complex<double> CubicFormula : : cubicRoot1 ( ) {
_u=_uu ; _v=_vv ;
r e turn ( _u+_v ) ;
}
complex<double> CubicFormula : : cubicRoot2 ( ) {
_u=_uu*CROOT2 ;
_v=_vv*CROOT3 ;
r e turn ( _u+_v ) ;
}
complex<double> CubicFormula : : cubicRoot3 ( ) {
_u=_uu*CROOT3 ;
_v=_vv*CROOT2 ;
r e turn ( _u+_v ) ;
}
/*Full cubicFormula **************************************/
186
FullCubicFormula : : FullCubicFormula ( ) {
// default , Pascal ' s triangle , a0 xˆ3+a xˆ2+b x+c==0;
CROOT2=cartesian (−0.5 ,0 .8660254037844386) ;//−1/2+i* s q r t (3) /2
CROOT3=cartesian (−0.5 ,−0.8660254037844386) ;//−1/2−i* s q r t (3) /2
setFormula ( 1 , 3 , 3 , 1 ) ; // solution x=−1 ( triple )
}
FullCubicFormula : : FullCubicFormula ( double a0 , double a , double b , double c ){
//a0 xˆ3+a xˆ2+b x+c==0;
CROOT2=cartesian (−0.5 ,0 .8660254037844386) ;//−1/2+i* s q r t (3) /2
CROOT3=cartesian (−0.5 ,−0.8660254037844386) ;//−1/2−i* s q r t (3) /2
setFormula ( a0 , a , b , c ) ;
}
int FullCubicFormula : : setFormula ( double a0 , double a , double b , double c ){
int returnValue =0;
isCubic_=1;
i f ( a0==0){
isCubic_=0;// f l a g
// cout<<”Warning : Cubic term vanishes !”<<endl ;
FullQuadraticFormula fullQuadraticTemp (a , b , c ) ;
// degenerate to xˆ2 formular
root1_=fullQuadraticTemp . quadraticRoot1 ( ) ;
root2_=fullQuadraticTemp . quadraticRoot2 ( ) ;
root3_=root2_ ;
returnValue=−1;
}
e l s e {
a/=a0 ; // normalize
b/=a0 ;
c/=a0 ;
a_=a /3 . 0 ; // convert to canonic form
p_=b−a*a / 3 . 0 ;
187
q_=c+(2.0* cube (a ) −9.0*a*b ) / 2 7 . 0 ;
discriminant_=cube ( p_ /3 . 0 )+square ( q_ /2 . 0 ) ;
// http : // en . wikipedia . org/wiki/Cubic_equation
i f ( discriminant_==0){
sqrtD_=0;
uu_=cbrt(−q_ /2 . 0 ) ;
vv_=cbrt(−q_ /2 . 0 ) ;
returnValue =0;
}
i f ( discriminant_ >0){
sqrtD_=sqr t ( discriminant_ ) ;
uu_=cbrt ( r e a l (−q_/2.0+ sqrtD_ ) ) ;
vv_=cbrt ( r e a l (−q_/2.0− sqrtD_ ) ) ;
returnValue =1;
}
e l s e {
sqrtD_=cartesian (0 , s q r t (−discriminant_ ) ) ;
uu_=cbrt(−q_/2.0+ sqrtD_ ) ;
vv_=cbrt(−q_/2.0− sqrtD_ ) ;
returnValue =3;
}
r e turn returnValue ;
}
}
int FullCubicFormula : : realRootNumber ( ) {
i f ( isCubic_ ){
i f ( discriminant_ <0)
r eturn 3 ;
e l s e
r e turn 1 ;
}
188
e l s e
r e turn −1;
}
complex<double> FullCubicFormula : : cubicRoot1 ( ) {
i f ( isCubic_ ){
u_=uu_ ; v_=vv_ ;
r e turn ( u_+v_ )−a_ ;
}
e l s e
r e turn root1_ ;
}
complex<double> FullCubicFormula : : cubicRoot2 ( ) {
i f ( isCubic_ ){
u_=uu_*CROOT2 ;
v_=vv_*CROOT3 ;
r e turn ( u_+v_ )−a_ ;
}
e l s e
r e turn root2_ ;
}
complex<double> FullCubicFormula : : cubicRoot3 ( ) {
i f ( isCubic_ ){
u_=uu_*CROOT3 ;
v_=vv_*CROOT2 ;
r e turn ( u_+v_ )−a_ ;
}
e l s e
r e turn root3_ ;
}
189
FDTDtime.h
#ifndef FDTDtime_h
#define FDTDtime_h
/* Estimate finishing time of the simulation program */
int FDTDtimer ( int total_I , int current_i , int inc_i , clock_t previous_clock ) ;
#endif
190
FDTDtime.cpp
#include ”FDTDarray . h”
#include <ctime>
#include <iostream>
using namespace std ;
/*−−−−−−−−−−−−−−−−−−−−Timer and/or estimated time f unct i on−−−−−−−−−−−−*/
int FDTDtimer ( int total_I , int current_i , int inc_i , clock_t previous_clock ) {
int i=current_i ;
int totalII=total_I ;
time_t ptimer2 ;
time_t ptimer3 ;
clock_t clock1=previous_clock ;
clock_t clock2=clock ( ) ; // current c l ock
clock_t clock3 ;
clock_t clock_diff ;
clock_t clock_total ;
struct tm *timeinfo2 ;
struct tm *timeinfo3 ;
int second_total ;
int second_frac ;
int hour_total ;
int minute_total ;
cout<<endl ;
//cout<<”i=”<<i<<endl ;
191
time(&ptimer2 ) ;
//cout<<”ptimer2===”<<ptimer2<<endl ;
//cout<<”clock1===”<<clock1<<endl ;
//cout<<”clock2===”<<clock2<<endl ;
clock_diff=clock2−clock1 ;
clock_total=(int ) ( totalII−(double )i ) /inc_i*clock_diff ;
//cout<<”clock_total tick remaining . . . ( clock_t )”<<clock_total<<endl ;
second_total=clock_total /CLOCKS_PER_SEC ;
second_frac=second_total ;
//cout<<”second_total remaining . . . ( s ec )”<<second_total<<endl ;
minute_total=second_frac /60 ;
second_frac−=minute_total *60 ;
hour_total=minute_total /60 ;
minute_total−=hour_total *60 ;
cout<<”Total time estimated : ”<<hour_total<<” hours , ”<<
minute_total<<” mintes and ”<<second_frac<<” seconds˜˜”<<endl ;
//cout<<”c l ock loop d i f f ( s ec )=”<<(double ) clock_diff /CLOCKS_PER_SEC<<endl ;
//cout<<”current time ( time_t )=”<<ptimer2<<endl ;
//cout<<endl ;
clock3=clock2+clock_total ;
//cout<<”c l ock now ( clock_t )=”<<c l ock ( )<<endl ;
//cout<<”finishing time ( clock_t )”<<clock3<<endl ;
//cout<<endl ;
ptimer3=second_total+ptimer2 ;
//cout<<”finishing timer ( time_t )===”<<ptimer3<<endl ;
192
timeinfo2=localtime (&ptimer2 ) ;
//cout<<”current asctime===”<<asctime ( timeinfo2 )<<endl ;
timeinfo3=localtime (&ptimer3 ) ;
cout<<”Finishing asctime :\ t”<<asctime ( timeinfo3 )<<endl ;
cout<<endl<<endl ;
r e turn 0 ;
}
193
REFERENCES CITED
[1] A. B. Kahng, Design Test of Computers, IEEE 27, 86 (2010), ISSN 0740-
7475.
[2] T.-Y. Chiang, B. Shieh, and K. C. Saraswat, IEEE Symp. VLSI Circuits, Dig.
Tech pp. 38–39 (2002).
[3] E. Ozbay, Science 311, 189 (2006), http://www.sciencemag.org/content/
311/5758/189.full.pdf, URL http://www.sciencemag.org/content/
311/5758/189.abstract.
[4] R. Zia, J. Schuller, A. Chandran, and M. Brongersma, Materials
Today 9, 20 (2006), ISSN 13697021, URL http://dx.doi.org/10.1016/
S1369-7021(06)71572-3.
[5] I. P. Kaminow, T. Li, and A. E. Willner, Optical fiber Telecommunications V
A: Components and Subsystems, Volume 1 (Academic Press, 2008).
[6] D. Pines, Rev. Mod. Phys. 28, 184 (1956).
[7] U. Fano, Phys. Rev. 103, 1202 (1956).
[8] S. Cunningham, A. Maradudin, and R. Wallis, Phys. Rev. B 10 (1974).
[9] P. Berini, R. Charbonneau, N. Lahoud, and G. Mattiussi, 98, 043109 (2005),
ISSN 00218979, URL http://dx.doi.org/doi/10.1063/1.2008385.
[10] W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature 424 (2003).
[11] M. L. Brongersma, J. W. Hartman, and H. A. Atwater, Phys. Rev. B
62, R16356 (2000), URL http://link.aps.org/doi/10.1103/PhysRevB.
62.R16356.
[12] K. Hasegawa, J. No¨ckel, and M. Deutsch, Appl. Phys. Lett. 84, 1835 (2004).
[13] F. Yang, J. R. Sambles, and G. W. Bradberry, Phys. Rev. Lett. 64, 559
(1990), URL http://link.aps.org/doi/10.1103/PhysRevLett.64.559.
194
[14] F. Yang, J. R. Sambles, and G. W. Bradberry, Phys. Rev. B 44, 5855 (1991),
URL http://link.aps.org/doi/10.1103/PhysRevB.44.5855.
[15] S. I. Bozhevolnyi and F. A. Pudonin, Phys. Rev. Lett. 78, 2823 (1997), URL
http://link.aps.org/doi/10.1103/PhysRevLett.78.2823.
[16] T. Nikolajsen, K. Leosson, and S. I. Bozhevolnyi, Applied Physics Letters 85
(2004).
[17] K. Hasegawa, J. No¨ckel, and M. Deutsch, Phys. Rev. A 75, 063816 (2008).
[18] M. Willatzen, Phys. Rev. A 80, 043805 (2009), URL http://link.aps.org/
doi/10.1103/PhysRevA.80.043805.
[19] S. A. Maier, M. L. Brongersma, P. G. Kik, S. Meltzer, A. A. G.
Requicha, and H. A. Atwater, Advanced Materials 13, 1501 (2001),
ISSN 1521-4095, URL http://dx.doi.org/10.1002/1521-4095(200110)
13:19<1501::AID-ADMA1501>3.0.CO;2-Z.
[20] J. A. Dionne, L. A. Sweatlock, H. A. Atwater, and A. Polman, Phys. Rev.
B 72, 075405 (2005), URL http://link.aps.org/doi/10.1103/PhysRevB.
72.075405.
[21] J. A. Dionne, K. Diest, L. A. Sweatlock, and H. A. Atwater, Nano Letters
9, 897 (2009), pMID: 19170558, http://pubs.acs.org/doi/pdf/10.1021/
nl803868k, URL http://pubs.acs.org/doi/abs/10.1021/nl803868k.
[22] K. F. MacDonald, Z. L. Sa´mson, M. I. Stockman, and N. I. Zheludev, Nature
Photonics 3, 55 (2009), URL http://www.nature.com/nphoton/journal/
v3/n1/pdf/nphoton.2008.249.pdf.
[23] N. Rotenberg, J. N. Caspers, and H. M. van Driel, Phys. Rev. B 80, 245420
(2009), URL http://link.aps.org/doi/10.1103/PhysRevB.80.245420.
[24] Z. L. Sa´mson, K. F. MacDonald, and N. I. Zheludev, Journal of Optics A:
Pure and Applied Optics 11, 114031 (2009), URL http://stacks.iop.org/
1464-4258/11/i=11/a=114031.
[25] M. Piliarik, L. Pa´rova´, and J. Homola, Biosensors and Bioelectronics 24
(2009).
[26] C. Nylander, B. Liedberg, and T. Lind, Sensors and Actuators 3, 79 (1982).
[27] B. Liedberg, C. Nylander, and I. Lundstro¨m, Sensors and Actuators 4, 299
(1983).
[28] J. Homola, S. S. Yee, and G. Gauglitz, Sensors and Actuators B 54, 3 (1999).
195
[29] B. Liedberg, C. Nylander, and I. Lundstro¨m, Biosensors Bioelectron 10, i
(1995).
[30] M. A. Cooper, Nature 1 (2002).
[31] R. Harris and J. Wilkinson, Sensors and Actuators B 29, 261 (1995).
[32] S. Ushioda and Y. Sasaki, Phys. Rev. B 27, 1401 (1983), URL http://link.
aps.org/doi/10.1103/PhysRevB.27.1401.
[33] A. Otto, I. Mrozek, H. Grabhorn, and W. Akemann, Journal of
Physics: Condensed Matter 4, 1143 (1992), URL http://stacks.iop.org/
0953-8984/4/i=5/a=001.
[34] A. C. R. Pipino, G. C. Schatz, and R. P. Van Duyne, Phys. Rev. B 49, 8320
(1994), URL http://link.aps.org/doi/10.1103/PhysRevB.49.8320.
[35] K. A. O’Donnell, R. Torre, and C. S. West, Phys. Rev. B 55, 7985 (1997),
URL http://link.aps.org/doi/10.1103/PhysRevB.55.7985.
[36] L. Cao, N. C. Panoiu, and R. M. Osgood, Phys. Rev. B 75, 205401 (2007),
URL http://link.aps.org/doi/10.1103/PhysRevB.75.205401.
[37] E. D. Palik, Handbook of Optical Constants of Solids (Academic Press, 1985).
[38] J. Jackson, Classical Electrodynamics (John Wiley & Sons, Inc., 1998).
[39] N. Bloembergen, Nonlinear Optics (World Scientific Publishing Co Pte Ltd,
1996).
[40] Y.-R. Shen, The Principles of Nonlinear Optics (John Wiley & Sons, 1984).
[41] R. W. Boyd, Nonlinear Optics (Academic Press, 2008).
[42] Y. Guo and M. Deutsch, Frontiers in Optics (FiO)/Laser Science XXVI (LS)
conference 696 (2010).
[43] Y. Guo and M. Deutsch, APS March Meeting (2010).
[44] S. Maier and H. A. Atwater, Journal of applied physics 98 (2005).
[45] P. Berini, Optics Letters 24 (1999).
[46] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Brooks Cole, 1976).
[47] S. Maier, Plasmonics - Fundamentals and Applications (Springer, 2007).
196
[48] M. J. Keskinen, P. Loschialpo, D. Forester, and J. Schelleng, 88, 5785 (2000),
ISSN 00218979, URL http://dx.doi.org/10.1063/1.1289045.
[49] H. S. Nalwa, Handbook of Advanced Electronic and Photonic Materials and
Devices: Conducting polymers (Elsevier, 2001).
[50] K. Lee, R. Menon, C. O. Yoon, and A. J. Heeger, Phys. Rev. B 52, 4779
(1995), URL http://link.aps.org/doi/10.1103/PhysRevB.52.4779.
[51] K. Lee and A. J. Heeger, Phys. Rev. B 68, 035201 (2003), URL http://
link.aps.org/doi/10.1103/PhysRevB.68.035201.
[52] S. R. Nagel and S. E. Schnatterly, Phys. Rev. B 9, 1299 (1974), URL http:
//link.aps.org/doi/10.1103/PhysRevB.9.1299.
[53] P. G. Etchegoin, E. C. L. Ru, and M. Meyer, 125, 164705 (2006), ISSN
00219606, URL http://dx.doi.org/doi/10.1063/1.2360270.
[54] I. R. Hooper and J. R. Sambles, Phys. Rev. B 65, 165432 (2002).
[55] D. Rogovin and T. P. Shen, Phys. Rev. B 37, 1121 (1988), URL http:
//link.aps.org/doi/10.1103/PhysRevB.37.1121.
[56] S. Palomba and L. Novotny, Phys. Rev. Lett. 101, 056802 (2008), URL http:
//link.aps.org/doi/10.1103/PhysRevLett.101.056802.
[57] C. Tang, Z. Wang, W. Zhang, S. Zhu, N. Ming, G. Sun, and P. Sheng,
Phys. Rev. B 80, 165401 (2009), URL http://link.aps.org/doi/10.1103/
PhysRevB.80.165401.
[58] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, Phys. Rev. Lett. 7,
118 (1961), URL http://link.aps.org/doi/10.1103/PhysRevLett.7.118.
[59] F. Moulin and D. Bernard, Optics Communications 164, 137 (1999), ISSN
0030-4018, URL http://www.sciencedirect.com/science/article/pii/
S0030401899001698.
[60] D. F. Eaton, Science 253, 281 (1991).
[61] M. J. Comstock, Materials for Nonlinear Optics (American Chemical Society,
19911).
[62] J. Zyss,Molecular nonlinear optics: Materials, physics and devices (Academic
Press, 1997).
[63] H. S. Nalwa and S. Miyata, Nonlinear optics of organic molecules and
polymers (CRC Press, 1997).
197
[64] N. Bloembergen, Applied Physics B 68 (1999).
[65] J. E. Sipe, V. C. Y. So, M. Fukui, and G. I. Stegeman, Phys. Rev. B 21, 4389
(1980), URL http://link.aps.org/doi/10.1103/PhysRevB.21.4389.
[66] F. De Martini, P. Ristori, E. Santamato, and A. C. A. Zammit, Phys. Rev. B
23, 3797 (1981), URL http://link.aps.org/doi/10.1103/PhysRevB.23.
3797.
[67] J. Renger, R. Quidant, N. van Hulst, S. Palomba, and L. Novotny, Phys.
Rev. Lett. 103, 266802 (2009), URL http://link.aps.org/doi/10.1103/
PhysRevLett.103.266802.
[68] D. Moss, J. E. Sipe, and H. M. van Driel, Physical Review B 35, 1130 (1986).
[69] G. Lu¨pke and G. Marowsky, Applied Physical B 53, 71 (1991).
[70] H. J. Simon, D. E. Mitchell, and J. G. Watson, Phys. Rev. Lett. 33, 1531
(1974), URL http://link.aps.org/doi/10.1103/PhysRevLett.33.1531.
[71] C.-C. Tzeng and J. T. Lue, Phys. Rev. A 39, 191 (1989), URL http://link.
aps.org/doi/10.1103/PhysRevA.39.191.
[72] G. S. Agarwal and S. S. Jha, Phys. Rev. B 26, 482 (1982), URL http:
//link.aps.org/doi/10.1103/PhysRevB.26.482.
[73] M. G. Weber, Phys. Rev. B 33, 6775 (1986), URL http://link.aps.org/
doi/10.1103/PhysRevB.33.6775.
[74] J. Nkoma, R. Loudon, and D. R. Tilley, J. Phys. C: Solid State Phys. 7, 3547
(1974).
[75] H. E. J. Dodd, R.K.and Morris, Soliton and nonlinear wave equations
(Academic Press, 1982).
[76] N. J. Zabusky and M. D. Kruskal, Phys. Rev. Lett. 15, 240 (1965), URL
http://link.aps.org/doi/10.1103/PhysRevLett.15.240.
[77] C. S. Gardner, J. M. Greene, M. D. Kruskal, and R. M. Miura, Phys.
Rev. Lett. 19, 1095 (1967), URL http://link.aps.org/doi/10.1103/
PhysRevLett.19.1095.
[78] H. C. Yuen and B. M. Lake, 18, 956 (1975), ISSN 00319171, URL http:
//dx.doi.org/doi/10.1063/1.861268.
[79] M. Destrade and G. Saccomandi, Phys. Rev. E 73, 065604 (2006), URL
http://link.aps.org/doi/10.1103/PhysRevE.73.065604.
198
[80] S. Bednarek, B. Szafran, and K. Lis, Phys. Rev. B 72, 075319 (2005), URL
http://link.aps.org/doi/10.1103/PhysRevB.72.075319.
[81] W. Singhsomroje and H. J. Maris, Phys. Rev. B 69, 174303 (2004), URL
http://link.aps.org/doi/10.1103/PhysRevB.69.174303.
[82] I. V. Shadrivov, A. A. Sukhorukov, Y. S. Kivshar, A. A. Zharov, A. D.
Boardman, and P. Egan, Phys. Rev. E 69, 016617 (2004), URL http://
link.aps.org/doi/10.1103/PhysRevE.69.016617.
[83] A. R. Davoyan, I. V. Shadrivov, and Y. S. Kivshar, Opt. Express 17,
21732 (2009), URL http://www.opticsexpress.org/abstract.cfm?URI=
oe-17-24-21732.
[84] L. F. Mollenauer, R. H. Stolen, and J. P. Gordon, Phys. Rev. Lett. 45, 1095
(1980), URL http://link.aps.org/doi/10.1103/PhysRevLett.45.1095.
[85] A. Hasegawa and F. Tappert, Applied Physics Letters 23, 142 (1973), URL
http://link.aip.org/link/?APL/23/142/1.
[86] Y. S. Kivshar and G. P. Agrawal, Optical Solitons (Academic Press, 2003).
[87] R. H. Stolen and C. Lin, Phys. Rev. A 17, 1448 (1978), URL http://link.
aps.org/doi/10.1103/PhysRevA.17.1448.
[88] G. I. Stegeman and M. Segev, Science 286, 1518 (1999), URL http://www.
sciencemag.org/content/286/5444/1518.full.
[89] G. Stegeman, D. Christodoulides, and M. Segev, Selected Topics in Quantum
Electronics, IEEE Journal of 6, 1419 (2000), ISSN 1077-260X.
[90] D. Korteweg and G. de Vries, Philos. Mag. 39 (1895).
[91] J. R. Taylor, Optical Solitons - Theory and Experiment (Cambridge
University Press, 1992).
[92] F. Abdullaev, S. Darmanyan, and P. Khabibullaev, Optical Solitons
(Academic Press, 1982).
[93] A. D. Boardman, G. S. Cooper, A. A. Maradudin, and T. P. Shen, Physical
review B 34, 8273 (1986).
[94] G. Agrawal, Nonlinear Fiber Optics (Academic Press, 2006).
[95] M. Kyoung and M. Lee, Optics Communications 171, 145 (1999), ISSN
0030-4018, URL http://www.sciencedirect.com/science/article/pii/
S0030401899005519.
199
[96] R. del Coso and J. Solis, J. Opt. Soc. Am. B 21, 640 (2004), URL http:
//josab.osa.org/abstract.cfm?URI=josab-21-3-640.
[97] J.-H. Huang, R. Chang, P.-T. Leung, and D. P. Tsai, Optics Communications
282, 1412 (2009), ISSN 0030-4018, URL http://www.sciencedirect.com/
science/article/pii/S0030401808012480.
[98] D. Mihalache, G. I. Stegeman, C. T. Seaton, E. M. Wright, R. Zanoni, A. D.
Boardman, and T. Twardowski, Opt. Lett. 12, 187 (1987), URL http://ol.
osa.org/abstract.cfm?URI=ol-12-3-187.
[99] S. A. Ramakrishna and O. J. Martin, Opt. Lett. 30, 2626 (2005), URL http:
//ol.osa.org/abstract.cfm?URI=ol-30-19-2626.
[100] A. M. Weiner, J. P. Heritage, R. J. Hawkins, R. N. Thurston, E. M. Kirschner,
D. E. Leaird, and W. J. Tomlinson, Phys. Rev. Lett. 61, 2445 (1988), URL
http://link.aps.org/doi/10.1103/PhysRevLett.61.2445.
[101] Y. S. Kivshar and B. A. Malomed, Rev. Mod. Phys. 61, 763 (1989), URL
http://link.aps.org/doi/10.1103/RevModPhys.61.763.
[102] O. Katz, Y. Lahini, and Y. Silberberg, Optical Letters 33 (2008).
[103] K.-J. Boller, A. Imamolu, and S. E. Harris, Phys. Rev. Lett. 66, 2593 (1991),
URL http://link.aps.org/doi/10.1103/PhysRevLett.66.2593.
[104] A. Kumar, S. F. Yu, X. F. Li, and S. P. Lau, Opt. Express 16,
16113 (2008), URL http://www.opticsexpress.org/abstract.cfm?URI=
oe-16-20-16113.
[105] K. Hasegawa, C. Rohde, and M. Deutsch, Optics Letters 31 (2006).
[106] G. Veronis and S. Fan, J. Lightwave Technol. 25, 2511 (2007), URL http:
//jlt.osa.org/abstract.cfm?URI=jlt-25-9-2511.
[107] M. I. Stockman, Phys. Rev. Lett. 106, 156802 (2011), URL http://link.
aps.org/doi/10.1103/PhysRevLett.106.156802.
[108] E. Kreyszig, Advanced Engineering Mathematics (John Wiley & Sons, Inc,
1998).
[109] G. B. Arfken, H. J. Weber, and F. Harris, Mathematical Methods for
Physicists, Sixth Edition: A Comprehensive Guide (Academic Press, 2005).
[110] T. R. Taha and M. J. Ablowitz, Journal of Computational Physics 55, 203
(1984).
200
[111] C. Pozrikidis, Numerical Computation in Science and Engineering (Oxford
University Press, 1998).
[112] D.-C. ionescu and H. Igel, Lecture notes in computer science 2659, 807
(2003).
201