Copyright © 2009 American Scientific Publishers
All rights reserved
Printed in the United States of America
Journal of
Computational and Theoretical Nanoscience
Vol. 6, 1808–1826, 2009
Modeling of Plasmonic Waveguide
Components and Networks
Georgios Veronis1 ∗ , Şükrü Ekin Kocabaş2 , David A. B. Miller2 , and Shanhui Fan2
1
Department of Electrical and Computer Engineering and Center for Computation and Technology,
Louisiana State University, Baton Rouge, LA 70803, USA
2
Ginzton Laboratory, Stanford University, Stanford, CA 94305, USA
We review some of the recent advances in the simulation of plasmonic devices, drawing examples
from our own work in metal-insulator-metal (MIM) plasmonic waveguide components and networks.
We introduce the mode-matching technique for modeling of MIM waveguide devices. We derive
the complete set of orthogonal modes
that theby
MIM
waveguide
Delivered
Ingenta
to: supports and use it to apply the
mode-matching technique to the analysis
of
plasmonic
waveguide
networks. We also introduce
University of Houston
several different equivalent models for plasmonic waveguide components, such as the characteristic
IP : 129.7.158.43
impedance model for deep subwavelength MIM waveguides, the scattering matrix description of
Mon,
05 Oct 2009 14:12:14
MIM waveguide junctions, and equivalent circuit models. The model abstraction provided by these
equivalent models is important for the analysis and synthesis of device functions, as illustrated with
the design of a waveguide mode converter.
Keywords: Plasmonics, Surface Plasmons, Mode-Matching, Equivalent Circuits, Scattering
Matrix, Characteristic Impedance.
CONTENTS
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2. Challenges Involved in Modeling of Plasmonic Devices . .
3. General Purpose Simulation Methods for Plasmonic
Waveguide Devices . . . . . . . . . . . . . . . . . . . . . . . . .
3.1. Finite-Difference Frequency-Domain Method . . . . . .
3.2. Finite-Difference Time-Domain Method . . . . . . . . .
4. Mode-Matching Method . . . . . . . . . . . . . . . . . . . . . .
4.1. Spectrum of MIM Waveguides . . . . . . . . . . . . . . .
4.2. Mode Orthogonality . . . . . . . . . . . . . . . . . . . . .
4.3. Discretization of the Continuous Spectrum . . . . . . .
4.4. Mode Completeness . . . . . . . . . . . . . . . . . . . . .
4.5. Convergence of the Mode-Matching Method . . . . . .
5. Equivalent Models for Plasmonic Waveguide Components
5.1. Characteristic Impedance Model . . . . . . . . . . . . . .
5.2. Scattering Matrix Description of Junctions . . . . . . .
5.3. Circuit Model for the Waveguide Junction . . . . . . .
5.4. Cascade Connection of Junctions . . . . . . . . . . . . .
6. Summary and Conclusions . . . . . . . . . . . . . . . . . . . .
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . 1808
. . . 1810
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1812
1812
1813
1814
1814
1815
1817
1817
1817
1818
1819
1819
1821
1822
1824
1824
1. INTRODUCTION
Surface plasmons are electromagnetic waves that propagate along the interface of a metal and a dielectric. In
surface plasmons light interacts with the free electrons of
∗
Author to whom correspondence should be addressed.
1808
J. Comput. Theor. Nanosci. 2009, Vol. 6, No. 8
the metal which oscillate collectively in response to the
applied field.1 Recently, nanometer scale metallic devices
have shown the potential to manipulate light at subwavelength scales using surface plasmons. This could lead to
photonic circuits of nanoscale dimensions. The use of
nano-metallic structures could also bridge the size mismatch between modern electronic components with critical
dimensions on the order of tens of nanometers and the
micrometer scaled optical devices.2
Plasmonic waveguides have shown the potential to guide
subwavelength optical modes at metal–dielectric interfaces. Several different plasmonic waveguiding structures
have been proposed,3–8 such as metallic nanowires4 5 and
metallic nanoparticle arrays.6 7 Most of these structures
support a highly-confined mode only near the surface plasmon frequency. In this regime, the optical mode typically has low group velocity and short propagation length.
It has been shown however that a metal-insulator-metal
(MIM) structure with a dielectric region thickness of
∼100 nm supports a propagating mode with a nanoscale
modal size at a wavelength range extending from DC
to visible.9 Thus, such a waveguide could be potentially
important in providing an interface between conventional
optics and subwavelength electronic and optoelectronic
devices. There have been several theoretical studies of
MIM waveguides in the literature.9–18 Because of the
1546-1955/2009/6/1808/019
doi:10.1166/jctn.2009.1244
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
Georgios Veronis received the B.S. degree in electrical engineering from the National Technical University of Athens, Greece, in 1997, and the M.S. and Ph.D. degrees in electrical
engineering from Stanford University, Stanford, CA, in 1999 and 2002 respectively. He was
an engineering research associate at the Ginzton Laboratory at Stanford University. He is
currently an assistant professor jointly at the Department of Electrical & Computer Engineering and the Center for Computation & Technology at Louisiana State University, Baton
Rouge, LA. He has published more than 25 refereed journal papers, and is the holder of
1 U.S. patent. His research interests include the theoretical analysis of nanophotonic and
plasmonic devices, and computational electromagnetics. Dr. Veronis is a member of IEEE
and OSA.
Şükrü Ekin Kocabaş received the B.S. and M.S. degrees in electrical engineering from
Bilkent University, Ankara, Turkey and Stanford University, Stanford, CA, USA in 2002
and 2004 respectively. Currently, he is a Ph.D. student majoring in electrical engineering
with a minor in physics. He has been a member of IEEE since 2000. His main research
interest is modeling and fabrication of nanophotonic optoelectronic devices.
Delivered by Ingenta to:
University of Houston
IP : 129.7.158.43
Mon, 05 Oct 2009 14:12:14
David A. B. Miller received the B.Sc. degree from St Andrews University, and the Ph.D.
degree from Heriot-Watt University, in 1979, both in Physics. He was with Bell Laboratories, from 1981 to 1996, as a department head from 1987, latterly of the Advanced Photonics
Research Department. He is currently the W. M. Keck Professor of Electrical Engineering at
Stanford University, Stanford, CA, the Director of the Solid State and Photonics Laboratory
at Stanford and a Co-Director of the Stanford Photonics Research Center. He also served
as Director of the Ginzton Laboratory at Stanford from 1997–2006. His research interests
include nanophotonic and quantum-confined optoelectronic physics and devices, and fundamentals and applications of optics in information sensing, switching, and processing. He
has published more than 200 scientific papers, delivered more than 100 conference invited
talks, and holds 62 patents. Dr. Miller has served as a Board member for both the Optical
Society of America (OSA) and the IEEE Lasers and Electro-Optics Society (LEOS), and in various other society and
conference committees. He was President of the IEEE Lasers and Electro-Optics Society in 1995. He has also served on
Boards for various photonics companies. He was awarded the Adolph Lomb Medal and the R. W. Wood Prize from the
OSA, the International Prize in Optics from the International Commission for Optics, and the IEEE Third Millennium
Medal. He is a Fellow of the Royal Societies of London and Edinburgh, OSA, APS, and IEEE, and holds honorary
degrees from the Vrije Universiteit Brussel and Heriot-Watt University.
Shanhui Fan is an associate professor of electrical engineering at Stanford University. He
received his Ph.D. in 1997 in theoretical condensed matter physics from the Massachusetts
Institute of Technology (MIT), and was a research scientist at the Research Laboratory
of Electronics at MIT prior to his appointment at Stanford. His research interests are in
computational and theoretical studies of solid state and photonic structures and devices,
especially photonic crystals, microcavities, and nanophotonic circuits and elements. He has
published over 140 refereed journal articles, has given over 100 invited talks, and was
granted 28 US patents. Professor Fan received a National Science Foundation Career Award
(2002), a David and Lucile Packard Fellowship in Science and Engineering (2003), the
National Academy of Sciences Award for Initiative in Research (2007), and the Adolph
Lomb medal from the Optical Society of America (2007). Dr. Fan is fellow of OSA and
APS, a senior member of IEEE, and a member of SPIE.
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
1809
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
predicted attractive properties of MIM waveguides, people
2. CHALLENGES INVOLVED IN MODELING
have also started to explore such structures experimentally.
OF PLASMONIC DEVICES
In particular, Dionne et al.19 have recently demonstrated
Surface plasmons can be described by macroscopic
waveguiding in a quasi-two-dimensional MIM geometry
electromagnetic theory, i.e., Maxwell’s equations, if the
experimentally, showing clear evidence of a subwavelength
electron mean free path in the metal is much shorter than
guided mode with substantial propagation distances.
the plasmon wavelength. This condition is usually fulfilled
Most of the theoretical studies of nanoscale plasmonic
at optical frequencies.24 Assuming an expit harmonic
devices involve the use of general purpose electromagtime dependence of all field quantities, Maxwell’s curl
netic simulation techniques. These techniques are able
equations in the frequency domain take the form
to simulate devices of arbitrary geometry and material
composition. On the other hand, as described in the
× Er = −i0 Hr
(1)
next section, these techniques typically require numerical
× Hr = i 0 r rEr
(2)
grids with resolution of tens to thousands of grid points
per wavelength, depending on the application. They are
In macroscopic electromagnetic theory, bulk material proptherefore not suitable for the design and optimization of
erties, such as the relative dielectric constant r r, are
multi-component optical circuits. One approach to address
used to describe objects irrespective of their size. However,
this problem is the development of special-purpose simulafor particles of nanoscale dimensions a more fundamental
tion techniques with superior computational efficiency for
description of their optical and electronic properties may
a specific class of problems. Previously, for example,
in the by Ingenta to:25
Delivered
be required.
modeling of photonic crystal devices, it has been
shown of Houston
University
Analytical methods, such as Mie theory,26 can only be
that the use of Wannier functions results in severalIPorders
: 129.7.158.43
applied to planar geometries or to objects with very high
of magnitude speedup in the optimizationMon,
of photonic
05 Oct 2009
14:12:14
symmetry
(spheres, infinite cylinders) and have therefore
20
crystal circuits. A second approach is the development
limited importance in the analysis of plasmonic devices
of equivalent models with only a few dynamic variables,
and structures. Thus, the analysis of plasmonic devices is
which are nevertheless capable of describing complex optimostly based on numerical simulation techniques.
cal processes in photonic devices in detail. For example,
Numerical modeling of plasmonic devices involves sevin the modeling of photonic crystal devices, it has been
eral challenges specific to plasmonics which need to be
shown that such models can be derived with the use of
addressed. The dielectric constant of metals at optical
coupled mode theory,21 22 and perturbation theory.22 23
wavelengths is complex, i.e., r = Re + i Im
In this context, here, we provide a review of some of
and is a complicated function of frequency.27 Thus, sevour own recent research activities aiming to advance the
eral simulation techniques which are limited to lossless,
theory and simulation of plasmonic waveguide devices
non-dispersive materials are not applicable to plasmonic
through the development of efficient special-purpose simdevices. In addition, in time-domain methods the disperulation techniques and equivalent models. We first give
sion properties of metals have to be approximated by
an overview of the challenges involved in modeling of
suitable analytical expressions.28 In most cases the Drude
plasmonic devices, and we briefly examine two general
model is invoked to characterize the frequency dependence
purpose simulation techniques which are widely used for
of the metallic dielectric function29
modeling of plasmonic waveguide devices. We then intro2p
duce the mode-matching technique for modeling of MIM
(3)
r Drude = 1 −
+ i
waveguide devices. We derive the complete set of orthogonal modes that the MIM waveguide supports and use
where p , are frequency-independent parameters. Howit to apply the mode-matching technique to the analyever, the Drude model approximation is valid over a limsis of plasmonic waveguide networks. We show that this
ited wavelength range.29 The range of validity of the Drude
special purpose simulation technique is far more effimodel can be extended by adding Lorentzian terms to
cient for this class of problems than general purpose elecEq. (3) to obtain the Lorentz–Drude model29
tromagnetic simulation techniques. Finally, we introduce
k
several different equivalent models for plasmonic waveg
fj 2j
(4)
=
+
r
LD
r
Drude
2
uide components, such as the characteristic impedance
2
j=1 j − − i j
model for deep subwavelength MIM waveguides, the
scattering matrix description of MIM waveguide juncwhere j and j stand for the oscillator resonant
tions, and equivalent circuit models. We show that the
frequencies and bandwidths respectively, and fj are
model abstraction provided by these equivalent models is
weighting factors. Physically, the Drude and Lorentzian
important for the analysis and synthesis of device functerms are related to intraband (free-electron) and interband
tions, and we illustrate this with the design of a mode
(bound-electron) transitions respectively.29 Even though
converter.
the Lorentz–Drude model extends the range of validity of
1810
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
grid resolution is required at the metal-dielectric interanalytical approximations to metallic dielectric constants,
face to adequately resolve the local fields. In addition,
it is not suitable for description of sharp absorption edges
several plasmonic devices are based on components of
observed in some metals, unless a very large number of
subwavelength dimensions.1 In fact, most of the potential
terms is used.29 In particular, the Lorenz-Drude model
applications of surface plasmons are related to subwavecannot approximate well the onset of interband absorplength optics. The nanoscale feature sizes of plasmonic
tion in noble metals (Ag, Au, Cu) even if five Lorentzian
devices pose an extra challenge to numerical simulation
terms are used.29 In Figure 1 we compare the Drude and
techniques.
Lorentz–Drude models with experimental data for silver.
We illustrate the challenges involved in modeling plasWe observe that even a five-term Lorentz–Drude model
monic devices using a simple example. We consider an
with optimal parameters results in a factor of two error
infinite periodic array of silver cylinders illuminated by a
at certain frequencies. An alternative approach is the use
plane wave at normal incidence (inset of Fig. 2(a)). We use
of an analytical expression based on multiple complexthe finite-difference frequency-domain method, described
conjugate pole-residue pairs.30 It has been shown that,
in more detail below, to calculate the transmission of
when such an approach is used with time-domain meththe periodic array. This method allows us to directly use
ods, it can lead to substantial savings in both memory and
experimental data for the frequency-dependent dielectric
computation time.30
constant of metals, including both the real and imagiIn addition, in surface plasmons propagating along
nary parts, with no further approximation. The fields are
the interface of a metal and a dielectric, the field is
discretized
concentrated at the interface, and decays exponentially
Delivered by Ingenta
to:on a uniform two-dimensional grid with grid
size
x
=
away from the interface in both the metal and
dielecUniversity of Houston y = l. In Figure 2(a) we show the calculated transmission as a function of frequency. We also
tric regions.1 Thus, for numerical methods based IP
on: dis129.7.158.43
show
the transmission of the structure calculated with
cretization of the fields on a numerical grid,
a
very
fine
Mon, 05 Oct 2009
14:12:14
(a)
(a)
0.9
0.8
102
k0
0.7
|εRe|
Transmission
2a
101
|εIm|
100
0.6
0.5
a
0.4
0.3
0.2
Drude
0.1
10−1
0
200
400
600
800
200
300
1000
400
500
600
700
800
900
Frequency (THz)
Frequency (THz)
(b)
(b)
0.85
102
|εRe|
Transmission
0.8
1
10
|εIm|
100
0.75
0.7
0.65
0.6
−1
10
0.55
0.5
200
400
600
800
1000
1
2
5
10
∆l (nm)
Frequency (THz)
Fig. 1. Real and imaginary part of the dielectric constant of silver at
optical frequencies. The solid lines show experimental data.27 The dashed
lines show values calculated using (a) the Drude model, (b) the Lorentz–
Drude model with five Lorentzian terms. The parameters of the models
are obtained through an optimization procedure.29
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Fig. 2. (a) Calculated transmission spectrum of an infinite array of silver cylinders (shown in the inset) for normal incidence and TM polarization. Results are shown for a = 100 nm. The dashed line shows the
transmission spectrum calculated using the Drude model with parameters
p = 137 × 1016 sec−1 , = 729 × 1013 sec−1 . (b) Calculated transmission
at 855 THz as a function of the spatial grid size l.
1811
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
the Drude model of Eq. (3). We observe that the use
with node coordinates rijk = xi yj zk is the simplest and
most commonly-used. A field quantity at nodal location
of the Drude model results in substantial error. In genrijk is denoted for convenience as fijk = f rijk . Based on
eral, the Drude model parameters are chosen to minimize
Eq. (5), the first derivative can be approximated by the
the error in the dielectric function in a given frequency
following central-difference formula
range.31 However, this approach gives accurate results in
a limited wavelength range, as illustrated in this example.
df fi+1 − fi−1
In general, the complicated dispersion properties of met≃
(6)
dx i
2x
als at optical frequencies pose a challenge in modeling of
plasmonic devices not encountered in modeling of low- or
which is second-order accurate, based on the discussion
high-index-contrast dielectric devices.
above, if the rectangular grid is uniform, i.e., xi = ix.
In Figure 2(b) we show the calculated transmission at
Similarly, the second derivative can be approximated by
a specific wavelength of 0 = 351 nm as a function of
the formula
the spatial grid size l. We observe that a grid size of
d 2 f fi+1 − 2fi + fi−1
l ≃ 1 nm is required in this case to yield reasonably accu≃
(7)
dx2 i
x2
rate results. The required grid size is directly related to
the decay length of the fields at the metal-dielectric interwhich is also second-order accurate on a uniform grid.28
face. In general, modeling of plasmonic devices requires
By replacing derivatives in differential equations with
much finer grid resolution than modeling of low- or hightheir finite-difference approximations, we obtain algebraic
index-contrast dielectric devices, due to the high
localiza- by Ingenta
Delivered
equationsto:
which relate the value of the field at a spetion of the field at metal-dielectric interfaces ofUniversity
plasmonic ofcific
node to the values at neighboring nodes. To solve
Houston
devices. The required grid size depends on the shape
Maxwell’s equations with the FDFD method, we discretize
IP : and
129.7.158.43
feature size of the modeled plasmonic device,
the
the 14:12:14
system of the three coupled scalar partial differential
Mon,metallic
05 Oct 2009
material used and the operating frequency.
equations obtained from the wave equation for the electric
field
3. GENERAL PURPOSE SIMULATION
METHODS FOR PLASMONIC
WAVEGUIDE DEVICES
Here we briefly examine two general purpose electromagnetic simulation techniques which are widely used for modeling of plasmonic waveguide devices, the finite-difference
frequency-domain (FDFD) and the finite-difference timedomain (FDTD) methods. We also examine how they
address the challenges mentioned above. Other general
purpose simulation techniques for modeling of plasmonic
waveguide devices include the Green dyadic method,4 the
finite-element method,32 and the method of lines.33
3.1. Finite-Difference Frequency-Domain Method
In finite-difference methods, derivatives in differential
equations are approximated by finite differences. To
approximate the derivative df /dxx0 we consider Taylor
series expansions of f x about the point x0 at the points
x0 + x and x0 − x and obtain28
df
f x0 + x − f x0 − x
=
+ Ox2 (5)
dx x0
2x
Equation (5) shows that a central-difference approximation
of the first derivative is second-order accurate, meaning
that the remainder term in Eq. (5) approaches zero as the
square of x.
In finite-difference methods a continuous problem is
approximated by a discrete one. Field quantities are
defined on a discrete grid of nodes. The rectangular grid
1812
× × Er − r r
2
Er = −i0 Jr
c2
(8)
For simplicity we consider here two-dimensional problems
with TE polarization. For TE polarization we have E = Ez ẑ
and the wave equation for the electric field becomes32 34
2
2
2
+
+
k
x
y
Ez x y = −i0 Jz x y (9)
0 r
x2 y 2
For simplicity we consider a uniform rectangular grid with
xi = ix, yj = jy, and replace the derivatives in Eq. (9)
with their finite-difference approximations of Eq. (7) to
obtain
fi+1j −2fij +fi−1j fij+1 −2fij +fij−1
+
+k02 rij fij = Aij
x2
y2
(10)
where f = Ez and A = −i0 Jz . Thus, application of
finite-difference approximations at the node location rij =
xi yj results in a linear algebraic equation which relates
the field fi j to the fields at the four adjacent nodes
fi+1 j , fi−1 j , fi j+1 , fi j−1 . By applying the finite-difference
approximation to all nodes of the grid we obtain a system of linear equations of the form Ax = b, where b is
determined by the source current J. Since the equation for
the field at each point involves only the fields at the four
(six in three dimensions, two in one dimension) adjacent
points, the resulting system matrix is extremely sparse.34
FDFD can be used to model plasmonic devices with
arbitrary geometries. In addition, FDFD is a frequencydomain technique and can thus treat arbitrary material
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
dispersion. Nonuniform and/or nonorthogonal grids are
required in FDFD for efficient treatment of curved surfaces and rapid field decays at metal-dielectric interfaces.
In FDFD, as in all other methods which are based on discretization of the differential form of Maxwell’s equations
in a finite volume, absorbing boundary conditions (ABCs)
are required so that waves are not artificially reflected at
the boundaries of the computational domain.28 32 Very efficient and accurate ABCs such as the perfectly matched
layer (PML) have been demonstrated.35 As mentioned
above, FDFD results in extremely sparse systems of linear equations. Such problems can be solved efficiently if
direct or iterative sparse matrix techniques are used.32 34
3.2. Finite-Difference Time-Domain Method
The FDTD method28 solves directly Maxwell’s timedependent curl equations
We observe that, using the Yee lattice, all spatial finitedifference expressions are central and therefore secondorder accurate. Similar finite-difference equations are
obtained by discretizing the other components of Eqs. (11)
and (12). In summary, FDTD is an explicit numerical
scheme which is second-order accurate both in time and
in space (in uniform media).
3.2.1. Treatment of Dispersive Media in FDTD
One of the major challenges in FDTD modeling of metals
at optical frequencies is the treatment of the metallic dispersion properties. As mentioned above, in time-domain
methods the dielectric constants of dispersive media have
to be approximated by suitable analytical expressions. The
most common algorithm for modeling dispersive materials
with FDTD is the auxiliary differential equation (ADE)
method.28 36 In dispersive materials relates E and D
H
Delivered
D = E
(15)
(11) by Ingenta to:
t
University of Houston
ADE is based on integrating an ordinary differential equaIP : 129.7.158.43
E
×H = 0 r
(12)
tion14:12:14
in time that relates Dt to Et, concurrently with
Mon, 05 Oct 2009
t
Maxwell’s equations. This equation is derived by taking
so that both space and time have to be discretized.
the inverse Fourier transform of Eq. (15).
The standard FDTD is based on the Yee algorithm.28
We consider here a simple example where the dielectric
As mentioned in the previous section, central-difference
constant r consists of a single Lorentzian term, i.e.,
approximations are second-order accurate. To achieve
second-order accuracy in time, the Yee algorithm uses a
20
=
(16)
28
r
leapfrog arrangement. E fields are calculated at t = nt
20 − 2 − i 0
using previously calculated and stored H fields. Then H
If we substitute Eq. (16) into Eq. (15) and take the inverse
fields are calculated at t = n + 1/2t using the preFourier transform, we obtain a second-order differential
viously calculated and stored E fields, and the process
equation relating D and E
continues until time-stepping is concluded. Applying this
scheme to Eq. (11) we obtain
D 2 D
+ 2 = 20 0 E
(17)
20 D + 0
t
t
t
n+1/2
n−1/2
n
H
= H
− × E
(13)
0
Equation (17) is discretized using a second-order accurate
central-difference scheme similar to those described above.
We observe that the leapfrog scheme yields centralWe note that, if the ADE method is used, E is obtained
difference in time and therefore second-order accurate
from
H in two steps. First, D is obtained from H by solvapproximations. In addition, since E(H) fields are obtained
ing
the
finite-difference approximation of
from previously calculated and stored H(E) fields, the
time-stepping is fully explicit, meaning that we do not
D
×H =
(18)
have to solve a system of simultaneous equations.28
t
To achieve second-order accuracy in space, FDTD
Second, E is obtained from D by solving the finiteuses a special grid, known as the Yee lattice, where
difference
approximation of Eq. (17). Calculation of the
every E component is surrounded by four H compofinite-difference
expressions of the first and second time
nents and every H component is surrounded by four E
28
derivatives
of
D
in Eq. (17) requires the storage of 2 precomponents. Based on this arrangement, discretization of
vious
values
of
D.
In other words not only Dn+1 but also
the x-component of Eq. (11) gives
n
n−1
D and D
are required to obtain E from D.
Hx n+1/2
Another approach to model dispersive materials with
i j k
FDTD is the recursive convolution (RC) method.28 36
t Ey ni j k+1/2 − Ey ni j k−1/2
n−1/2
= Hx i j k +
FDTD is a finite-difference method, so its performance
0
z
in
modeling plasmonic devices is similar to the perfor
Ez ni j+1/2 k − Ez ni j−1/2 k
mance of FDFD. However, there are some major differ−
(14)
ences. First, as already mentioned above, in time-domain
y
× E = −0
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
1813
Modeling of Plasmonic Waveguide Components and Networks
methods the dispersion properties of metals have to be
approximated by suitable analytical expressions which
introduce substantial error in broadband calculations. In
addition, the implementation of the ADE or RC methods
requires additional computational cost and extra memory
storage.28 36 On the other hand, in FDTD it is possible to
obtain the entire frequency response with a single simulation by exciting a broadband pulse and calculating the
Fourier transform of both the excitation and the response.28
Veronis et al.
εm
x
z
a
εi
a′
x=0
(a)
x=a+h
= a′ + h′
x
4. MODE-MATCHING METHOD
h
εm
z
h′
In this section, we introduce the mode-matching technique for MIM waveguide devices. The mode-matching
technique37 is commonly used in the microwave and optia
cal domains.38–41 We derive the complete set of orthogonal
εi
a′
modes that the MIM waveguide supports and use it to
x=0
apply the mode-matching technique to the calculation of
(b)
scattering at the junction between two guides Delivered
with differ- by Ingenta to:
3. (a) Geometry for the even modes of the MIM waveguide. The
University ofFig.
Houston
ent cross sections.
4.1. Spectrum of MIM Waveguides
x = 0 plane contains a fictitious perfect electric conductor (PEC) to simIP : 129.7.158.43
plify the problem when dealing only with even TM modes of the guide.
Mon, 05 Oct 2009
14:12:14
This fictitious MIM waveguide is equivalent to an actual guide with an
The mode-matching method is based on expanding the
fields in terms of the modes of the waveguides. In this
context here we derive the modal structure (spectrum) of
the MIM waveguide.42 We will specifically focus on the
even modes of the waveguide, for which the transverse
magnetic (TM) field component is an even function of the
transverse coordinate, x. The reason why we focus on even
modes is that we will be analyzing the scattering of the
fundamental, even mode of the MIM waveguide—which
is also a TM mode10 11 15 —off of a symmetric junction
with a different sized MIM waveguide. Due to the symmetry of the problem at hand, even modes will be sufficient.
Evenness of the function is achieved by putting a fictitious
perfect electric conductor (PEC) at the x = 0 plane of the
waveguide, which forces the tangential electric field Ez to
be an odd function, and the magnetic field Hy to be an
even function of x. In other words, the modes of this fictitious waveguide with the PEC at x = 0 are mathematically
the same as the even modes of the actual waveguide of
interest, and so we will work with this hypothetical waveguide. The geometry is as shown in Figure 3(a). m refers to
the permittivity of the metal region and i of the insulator
region. At infrared frequencies, m (Fig. 1) is a complex
number with a large, negative real part and a relatively
small imaginary part (the sign of which is determined by
the time convention used, being negative for an expit
time dependence).
We begin with Maxwell’s equations (Eqs. (1), (2)). The
MIM waveguide is a two dimensional structure which does
not have any variation in the y direction. Therefore, we can
eliminate all the derivatives with respect to y in Maxwell’s
equations. Furthermore, our study will be based on the TM
1814
insulator thickness of 2a. The inset shows the equivalent symmetric junction of two MIM waveguides. Dashed line in the inset is the plane of symmetry, which is where the fictitious PEC layer is introduced. (b) Geometry
for mode matching. x = a + h plane of Figure 3(a) is terminated by a PEC
which leads to a discretization of the continuous spectrum.
modes which only have the Hy , Ex and Ez field components. Also, the uniformity of the waveguide in the z direction leads to exp−ikz z as the space dependence in z by
using the separation of variables technique for differential equations (kz may, however, be a complex number).
After simplifying the curl equations in (1), (2), we have the
following relationships between the different field
components
d
E x
dx z
ikz Hy x = i xEx x
ixHy x = ikz Ex x +
(19)
d
H x = i xEz x
dx y
Using these equations we get the following differential
equation for Hy
d 1 d
+ 2 x x Hy = kz2 Hy (20)
x
dx x dx
and since Ez 0 = 0 by the PEC wall at x = 0, the boundary
condition for Hy under (19) becomes dHy x/dxx=0 = 0.
4.1.1. Point Spectrum
The dispersion equation that should be solved in order to
find the kz values for the modes of the MIM waveguide is
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
derived by satisfying the continuity of tangential electric
and magnetic fields at material boundaries and applying
the boundary conditions. We refer the reader to Refs. [43
and 9, 11, 44–46] for the details. The eigenvectors &n
and the dispersion equation for the corresponding eigen2
of (20) for the even TM modes of the MIM
values kzn
waveguide are
cosh'i n x 0 < x < a
&n x = H0 cosh'i n a
(21)
e−'m n x−a
a<x<
tanh'i n a = −
'm n /
'i n /
2
2
2
kz
n = 'm n +
m
m
(22)
Re'2m , < 0 and Im'2m , = 0. These conditions can be
written in terms of kz, by using (26) as
2
2
Rekz
, −
m
<0
and
2
2
Imkz
, −
m
=0
Note that when (22) holds true, we have - = −1 in (25)
which makes (24) and (21) equivalent.
4.2. Mode Orthogonality
Orthogonality and completeness are two very valuable
properties of modes, which make the mode matching technique possible. Here we use the pseudo-inner product,
··, defined as50
i
= '2i n + 2
i
f g =
(23)
f xgxdx
0
It can be shown that two different eigenfunctions, &1 x
where Re'm n > 0 so that &n x does not diverge and is
2
and
&2 x, corresponding to two different eigenvalues kz
1
integrable. Here n is a discrete index for the eigenvalues
−1
50 51
2
and kz 2 are
Delivered by Ingenta
to: pseudo-orthogonal with x weight
and the eigenfunctions.
University of Houston
−1 &1 &2 = 0
(27)
IP : 129.7.158.43
Mon, 05 Oct 2009
14:12:14
From
(19) it can be seen that −1 &1 is proportional to
4.1.2. Continuous Spectrum
In this section, we show that a continuous spectrum exists
in the MIM waveguide. The utility of the continuous spectrum will be evident in the mode matching analysis.
The condition of square integrability of the modes can
be replaced by the weaker condition of finiteness of the
modes in their domain of definition.47 For the MIM waveguide, this would imply a non-zero, yet finite electromagnetic field at infinity. These infinite-extent and, therefore,
infinite energy, continuum modes (which can be normalized through the use of the Dirac delta distributions48 49 ) are
integrated to realize any physically possible finite energy
field configuration. In this respect, such an approach is similar to the well-known Fourier transform methods, where
finite energy functions are expanded in terms of the infinite
energy exponentials.
Constraining fields to be finite, instead of zero, at infinity leads to the following field profile
cosh'i, x
0<x<a
cosh'
i, a
+, x = H0
(24)
cosh'm , x − a
+ - sinh'm , x − a a < x <
-=
'i, /
'm , /
i
tanh'i, a
(25)
m
2
kz,
= '2m , + 2
m
= '2i, + 2
i
(26)
which is calculated very similarly to the case of a dielectric slab.43 Here , is a continuous index for different
functions in the continuous spectrum. For finite +, , the
arguments inside the hyperbolic functions for x > a in
(24), 'm , , should be purely imaginary which implies that
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
the transverse electric field component Ex of the mode.
Therefore, the orthogonality condition can also be written
as
0
Ex1 xHy2 xdx =
A
E1 r × H2 r · dA = 0
which is the well known modal orthogonality condition
proved by the Lorentz reciprocity theorem,52 where A
denotes the cross section of the waveguide.
One can directly verify (27) by integration and using
'2m 1 − '2m 2 = '2i 1 − '2i 2 which is a result of (23). The
following orthogonality conditions between the elements
of the point (&n ) and the continuous (+, ) spectrum can
similarly be proved
−1
−1
&n +, = 0
for all n and ,
+ +, = 0
for
,=
In the following sections, we will be working with fields
at the junction of two different waveguides. For notational
abbreviation we will use the following convention
i
/L R2
i
/L R2
e/L R2 = Ex i
h/L R2 = Hy i
where /L R2 is used to denote the modes of the left and
right side of the junction, which leads to the following
orthogonality condition
i
j
e/L R2 h/L R2 = 3ij 4/L R2
(28)
where 3ij is the Kronecker delta function and 4 is the
overlap integral of the electric and magnetic transverse
fields.
1815
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
Table I. Adjectives.
fields have an exp−ikz z dependence. The argument principle method gives us the 'm value for the modes. By
Signifier
Signified
using (23) we get the kz2 value. We then calculate kz2 1/2
Leaky
Re'm < 0
and choose the root which satisfies Imkz < 0.
Proper
Re'm > 0
In Figure 4 the spectrum of an idealized lossless silverImproper
Re'm = 0
like MIM waveguide is shown on the plane of '2m for
Forward
Rekz > 0
Backward
Rekz < 0
m = −143497 which is the real part of the permittivity
of silver at a wavelength of 1550 nm.27 59 There are
four real modes for 2a = /4—TM0 , TM2 , TM4 , TM6 —
After the classification and analysis of the MIM wavegindexed according to the number of zero crossings in Hy .
uide modes, we will now visualize different parts of its
There is also an infinite number of complex modes, which
spectrum by finding the zeros of the respective disperare those with eight and more zero crossings in the insusion equations through the use of the argument principle
lator region. These modes have a 'm with a positive real
method.53–58 We will use the adjectives in Table I to further
part that is rather small compared to the imaginary part—
differentiate between the modes.
this can also be deduced from the scale of the imaginary
Leaky modes are not normalizable and are not part of the
axis of Figure 4. The continuous spectrum is illustrated
spectrum. Proper modes can be normalized by the usual
by the thick line which corresponds to Re'2m < 0 and
integration and they form the point spectrum. Improper
Im'2m = 0. This line is also the branch cut of the square
modes can be normalized by using the Dirac delta funcroot function that is used to get 'm from '2m . The field proDelivered
by
Ingenta
to:
tions, 3x. They form the continuous spectrum. Forward
files of the modes in the insulator region, as shown in the
University
of
Houston
modes have a positive phase velocity, whereas the backinsets of Figure 4, look quite similar to the field profile of
: 129.7.158.43
ward modes have a negative phase velocity. WeIPdecide
the even modes of a parallel plate waveguide with a plate
05 Oct
14:12:14
on the sign of Rekz based on Imkz : By Mon,
definition,
all 2009
separation
of 2a.
2
modes are propagating in the +z direction. Therefore, in
'm = 0 is the bifurcation point for the point specthe limit z → , the fields should go to zero. Such a
trum when m is purely real. For positive '2m , the point
spectrum has real modes, whereas for negative '2m , the
behavior is possible only if Imkz is negative, since the
TM8
TM6
0.4
TMC1
2
Im(km
)/k02
0.2
TM8
TMC2
TM4
TMC1
TM2
0.0
TM6TM4TM0
TMC2
TM2
–0.2
–0.4
–800
TM0
–600
–400
–200
0
200
2
Re(km
)/k02
Fig. 4. Spectrum of the MIM waveguide for m = −143497 and 2a = /4 where = 1550 nm is the wavelength of operation. There are four real
modes and an infinite number of complex modes, all denoted with the • symbol. The thick line denotes the continuous spectrum. Due to the fact
that m is real, complex modes come in complex conjugate pairs. Insets show the Hy mode shapes in the x direction for the discrete spectrum (TM0
through TM8 ) and the continuous spectrum (TMC1 and TMC2 )—solid lines in the insets are the real part of the mode, dashed lines are the imaginary
part. The locations of the drawn continuous modes are shown by the • symbol. Modes in the continuous spectrum are purely oscillatory in the x
direction. Complex modes have a small decay, which is not visually apparent in the inset for TM8 .
1816
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Veronis et al.
point spectrum splits into two branches that are complex conjugates of one another. '2m = 0 corresponds to
kz2 = m k02 which then implies kz = −i m k0 —bounded
modes should have Imkz < 0.
When loss is introduced to the metal, the spectrum
moves on the complex plane.42 The forward, proper, complex modes of the lossless case turn into leaky modes by
migrating into the third quadrant of the complex 'm plane.
Modeling of Plasmonic Waveguide Components and Networks
waveguides convergence of the fields on both sides of the
junction is obtained only when the continuous spectrum
is also taken into consideration. Thus, it is clear that the
point spectrum on its own is not sufficient to describe the
behavior of the waveguide junctions. Inclusion of the continuous spectrum is essential. This is further illustrated in
the next section.
4.5. Convergence of the Mode-Matching Method
4.3. Discretization of the Continuous Spectrum
Now that we know how to treat the continuous spectrum
and are confident that the collection of the point and the
continuous spectrum results in a complete basis set, we
can proceed with the mode-matching formalism. We will
begin by assuming that the pth mode of the left waveguide propagates toward the right, scatters and creates the
following set of fields at the right and left sides of the
junction, which by the continuity of the tangential magnetic and electric fields, are set equal
The presence of a continuous spectrum leads to the formation of integral equations when the mode-matching method
is applied.48 The integral equation is then expanded
using an orthogonal basis set—not necessarily that of the
modes—to solve the scattering problem.
Another way to approach the scattering problem is
to limit the transverse coordinates by a PEC wall. This
approach has the effect of discretizing the continuum part
Delivered by Ingenta to:
of the spectrum47 60 —turning it into a discrete spectrum.
University of Houston
m
k
To limit parasitic reflections from the PEC walls, absorb3mp + Rmp hL x =
Tkp hR x
(31)
IP : 129.7.158.43
38
ing layers can be positioned before the PEC termination.
m=1
k=1
Mon, 05 Oct 2009 14:12:14
In Ref. 61 detailed analysis of how the continuous spec
m
k
3mp − Rmp eL x =
Tkp eR x
(32)
trum appears from a discrete collection can be found. We
m=1
k=1
will use a PEC wall to discretize the continuous spectrum.
Also, we will not use any perfectly matched layers to limit
Here Rmp is the reflection coefficient of the mth mode of
parasitic reflections since the metallic sections with perthe left waveguide in response to an incoming field in the
mittivity m effectively absorb the fields away from the
pth mode. Similarly, Tkp is the transmission coefficient of
junction.
the kth mode of the right waveguide. Note that we chose
The geometry is as shown in Figure 3(b). For the left
Rmp to denote the reflection coefficient for the transverse
waveguide the dispersion equation for modes becomes
magnetic fields, which automatically results in −Rmp as
the reflection coefficient for the transverse electric fields.
'm n / m
tanh'i n a = −
tanh'm n h
(29)
In Ref. [62], it is shown that the testing of the above
'i n / i
equations (i.e., the discretization of the equation using integration of both sides by a given function) should be done
which asymptotes to (22) as h → . The transverse magby the magnetic field of the larger waveguide for enforcing
netic field shape is
electric field continuity (32) and by the electric field of the
smaller waveguide to enforce the magnetic field continu cosh'i n x
0<x<a
cosh'i n a
ity (31). Although that analysis was specifically done for
&n x = H0
(30)
waveguides with perfect metals m → , we still use
cosh'm n x − a − h
that strategy so that the formulation limits to the correct
a < x < a+h
cosh'm n h
one should metals be made perfect.
For those cases where a < a′ , we will take the pseudon
n
inner product of (31) with eL and of (32) with hR . Fur4.4. Mode Completeness
thermore, assuming there are L modes on the left and R
modes on the right, we get
We tested the completeness of the modes by expanding the
R
L
fundamental mode of an MIM waveguide of a given thick
n
k
m
Tkp eL hR
3mp + Rmp 4L 3mn =
ness in terms of the modes of the MIM waveguide with a
m=1
k=1
different thickness. Without the continuous spectrum, the
R
L
field expansion converges, but to a field profile which is
k
m
n
=
Tkp 4R 3kn
3
−
R
e
h
42
mp
mp
R
L
not the same as the desired profile of the right junction.
m=1
k=1
On the other hand, inclusion of the continuous spectrum
through the discretization of the continuum by a PEC wall
with the help of (28). These are linear matrix equations
leads to the correct field profile. Similarly, we found that
with Rmp and Tkp as the unknowns. After calculating the
inner products, the set of equations can be inverted to
for the magnetic field profile at the junction of two MIM
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
1817
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
the mode-matching calculations to converge. Without the
continuous spectrum, the mode matching results converge
to the wrong result. Inclusion of the continuous spectrum
2a′
2a
decreases the error to around 2%, which is probably due
0.3
to the space discretization of FDFD simulations as well
as the method used in the de-embedding of the scattering
2a
=
0.9λ
0.2
coefficients from fields. As is also evident from Figure 5(a)
2a = 0.5λ
the utility of the single mode L = R = 1 mode-matching
0.1
calculations increases as the dimensions of the waveg2a = 0.1λ
uides decreases. The single mode approximation is closely
0.0
related to the characteristic impedance model described in
5
10
15
20
25
30
Number of modes
Section 5.1 which is a good approximation for deep subwavelength structures. In Figure 5(a) we also show the
90º
(b)
120º
effect of neglecting the backward modes in the mode60º
matching calculations for the 2a = 01 case. Backward
modes are important in this subwavelength geometry; however, for the wider geometries of the 2a = 05 and
150º
30º
2a = 09 cases we did not observe any increase in the
error when
Delivered by Ingenta
to:backward modes were neglected in the modematching
University of Houstoncalculations.
1
In Figure 5(b) we visualize the scattering coefficient
IP : 0º129.7.158.43
180º
S22
of
main mode of the MIM waveguide. (The extracMon, 05 Oct 2009the
14:12:14
tion of the scattering parameters from full field simula0.2
1
S11
tions is described in detail in the section below.) We do
0.4
the calculations in two different ways, one using FDFD,
330º
210º
and the other using the mode-matching technique with
0.6
the point and the continuous spectrum. When applying
0.8
mode-matching, we use the a > a′ formulation for S11
300º
240º
calculations and the a < a′ one for S22 . There is a very
1.0
good match between the results of the two techniques,
270º
verifying the applicability of the mode-matching method.
Fig. 5. (a) Convergence study of the reflection coefficient, S11 , of the
In addition, since the convergence of the mode-matching
main mode of the left waveguide traveling toward the right waveguide for
method requires relatively few modes, the use of the mode2a = 09, 2a = 05 and 2a = 01. a′ /a = 04 and m = −143497 −
i9517 for all cases. Dashed lines are for calculations including the point
matching technique results in orders of magnitude speedup
spectrum only. Solid lines are those with both the point and the conin the calculation of the properties of the junction.
(a)
Error
0.4
tinuous spectrum. Empty circles denote the calculations done with the
forward point spectrum and the continuous spectrum for the 2a = 01
MM
FDFD
FDFD
−S11
/S11
where MM stands for
case only. Error is defined as S11
mode-matching and FDFD for finite-difference frequency-domain calculations. The inset shows the junction geometry. (b) Reflection coefficient
of the main mode at the junction between two MIM waveguides of different insulator thicknesses, plotted on the complex plane within the unit
circle. Both waveguides have m = −143497−i9517 and i = 10. Filled
circles, •, are FDFD results, empty circles, , are mode-matching results.
The left waveguide thickness is fixed at 2a = 09. The right waveguide
thickness varies from 2a′ = /002 004 092. The origin is the
zero reflection point that corresponds to a = a′ . As a′ decreases progressively toward zero, we move progressively along the curves away from
the origin. The first set of curves, S11 , are for the case when the mode of
the left waveguide, traveling from left to right, is scattered by the junction. The second set of curves, S22 , are for the case when the main mode
of the right waveguide, traveling from right to left, is scattered by the
junction. Insets illustrate the respective cases.
give the reflection and transmission coefficients for the
modes.
In Figure 5(a), we compare the mode-matching method
with the FDFD technique. It takes relatively few modes for
1818
5. EQUIVALENT MODELS FOR PLASMONIC
WAVEGUIDE COMPONENTS
Modeling electromagnetic wave propagation using transmission lines has been one of the most important achievements of microwave network theory.63 The concept of
impedance64 and understanding the effects of waveguide
discontinuities in terms of lumped circuit elements were
crucial in this respect. Even though the properties of metals are quite different at optical wavelengths compared
to the microwave, designs that are qualitatively similar to
their low frequency counterparts have been demonstrated
at optical frequencies.65 It has also been shown that the
fundamental TM mode of an MIM waveguide continuously changes to the TEM mode of a parallel-plate waveguide with PEC boundaries as the frequency of operation
is decreased.66 It is intriguing to ask whether methods of
microwave can be applied to this new generation of nanometallic structures to come up with concise descriptions of
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Modeling of Plasmonic Waveguide Components and Networks
components that can lead to a simplified approach to the
design of functional systems composed of many interacting parts.
Here we first introduce the concept of the characteristic impedance for MIM plasmonic waveguides and we
show its validity and usefulness for subwavelength guides.
We then characterize the modal reflection and transmission
from MIM junctions using the scattering matrix approach,
a commonly used method in microwave network theory.
We also represent the scattering matrix of MIM junctions
in terms of an equivalent lumped circuit model. Finally, to
test our characterization, we design a cascade connection
of MIM junctions to couple the mode of a wavelengthsized MIM waveguide to that of a subwavelength one with
zero reflection. Throughout our analysis, we will compare
MIM waveguides to PEC parallel plate waveguides and
comment on the similarities and the differences between
the two.
0.12
0.1
0.08
Reflection
Veronis et al.
din
0.06
0.04
dout
0.02
0
1
1.5
2
2.5
3
din/dout
Fig. 6. Reflection coefficient R of a MIM T-shaped splitter (shown in
the inset) as a function of din /dout at 0 = 155 m calculated using
FDFD. We also show with dashed line the reflection coefficient R̄ calculated based on the characteristic impedance ZMIM and transmission-line
theory. Results are shown for dout = 50 nm.
Delivered by Ingenta to:
University ofcompared
Houstonto the wavelength and the quasistatic approxIP : 129.7.158.43
imation therefore breaks down. We found that similar
The characteristic impedance of the fundamental
Mon, 05TEM
Oct 2009
14:12:14
deviations
are observed for PEC parallel-plate waveguides.
5.1. Characteristic Impedance Model
mode in a PEC parallel-plate waveguide is uniquely
defined as the ratio of voltage V to surface current density
I and is equal to67
ZTEM ≡
V
E d :
= x = TEM d =
I
Hy
0
0
d
(33)
0
where Ex , Hy are the transverse components of the electric
and magnetic field respectively, and we assumed a unitlength waveguide in the y direction. For non-TEM modes,
such as the fundamental MIM mode, voltage and current
are not uniquely defined. However, metals like silver satisfy the condition metal ≫ diel at the optical communication wavelength of 1.55 m.27 Thus, Ex metal ≪ Ex diel
so that the integral of the electric field in the transverse
direction can be approximated by Ex diel d and we may
therefore define the characteristic impedance of the fundamental MIM mode as15
ZMIM d ≡
Ex diel d :MIM d
=
d
Hy diel
0
(34)
where :MIM d = 2;/g d, and g is the guide wavelength. In Figure 6 we show the reflection coefficient of
an MIM T-shaped splitter calculated based on ZMIM as
Z − Z0 2 2ZMIM dout − ZMIM din 2
=
R̄ = L
(35)
2Z d + Z d
ZL + Z0
MIM
out
MIM
in
We observe that there is very good agreement between
R̄ and the exact reflection coefficient R calculated using
FDFD. This agreement suggests that the concept of characteristic impedance for MIM waveguides is indeed valid
and useful. The deviation between R̄ and R at large values of din /dout is due to the fact that din is not very small
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Such deviations decrease at longer wavelengths in both the
PEC and MIM waveguide cases.
5.2. Scattering Matrix Description of Junctions
In this section, we characterize the modal reflection and
transmission from MIM junctions using the scattering
matrix approach, a commonly used method in microwave
network theory.68 We consider the geometry shown in
Figure 7(a). The insulating region is free space with a permittivity i = 1, and the metal is silver with a permittivity
of m = −143497 − i9517.27 59 The wavelength of operation is fixed at = 1550 nm, in the L band of optical
telecommunications.
Using the dispersion equation for even modes of the
MIM waveguide it can be shown that only a single even
propagating mode can exist for 2a < 097 for our choice
of m , i and . The condition for the PEC parallel plate
waveguide is similar, where only a single even propagation mode exists for 2a < 10. When there is only one
propagating mode, far away from the waveguide junction
the fields can be written in terms of that main mode of the
system since all higher order modes will have an exponential decay much faster compared to the main propagating
mode. Under such circumstances, the effects of the waveguide junction on the propagating modes can be described
using the single mode scattering matrix (S) formalism.69
In the terminology of the scattering matrix, the forward
and backward mode amplitudes are considered to scatter
from one “port" to another. Here we can think of the ports
as being the left and right port planes shown in Figure 7(a).
These ports are sufficiently far to the left and right of the
junction that the fields have settled down again to being
1819
Modeling of Plasmonic Waveguide Components and Networks
εm
5λ
5λ
Veronis et al.
1
1
S11
x
εi
2a
2a′
z
(b)
(c)
1
+
HL
Left port
Right port
1
S22
R
+
HR
S
–
(d)
–
HR
HL
(e)
(a)
Fig. 7. (a) Description of the modeling geometry. Dashed lines represent the location of the left and right ports of the overall scattering matrix S that
describes this junction (schematically shown in the bottom part of the figure). (b) Description of the steps taken in extracting S from fields. Calibration
simulations with uniform insulator widths of 2a and 2a′ , which give the wave vector k and the values of the incoming fields at the left, HL+ &L x,
Delivered by Ingenta to:
and right ports, HR+ &R x. (c) Field impinging from the left side, which leads to S11 . (d) Field impinging from the right side, which leads to S22 .
University
of Houston
(e) Simulation domain is terminated by a perfect electric conductor
at the right
input port plane. S12 is extracted from the reflection coefficient R using
IP : 129.7.158.43
the previously calculated S11 and S22 .
Mon, 05 Oct 2009 14:12:14
the propagating modes of the guides (nearer to the interface, there will in general be other field behavior, including various near-field components that decay rapidly with
distance).
If we can deduce the scattering matrix for such a junction, then we can have a very simple way of modeling
the behavior of structures containing such junctions, as is
already well known in the modeling of microwave guides.
The elements of the scattering matrix, S11 , S12 , S21 , S22 ,
are complex numbers which describe the phase and magnitude of the reflection and transmission of the main modes.
Thus, in general there are 8 independent real numbers
in S. However, under certain conditions the number of
independent parameters can be reduced. First of all, if the
system is composed of reciprocal media (i.e., symmetric permittivity and permeability tensors) then using the
Lorentz reciprocity theorem it can be shown that S12 = S21 .
Note that this equality implies a certain normalization of
the modes,69 specifically
A
EL × HL · dA =
A
ER × HR · dA = 1
(36)
where E/L R2 and H/L R2 denote the electric and magnetic
components of the main propagating modes on the left
(L) and the right (R) of the waveguide junction. Also note
that for lossless systems S is a unitary matrix69 (though
in general here we will be considering systems with loss).
As a result, using reciprocity it is possible to describe a
lossy junction using six real numbers, two for each of S11 ,
S12 , and S22 . When there is no loss we only need three real
numbers due to the unitarity of S.
We now describe the method we used to extract the elements of S from the electromagnetic fields in such waveguide junctions. We solved Maxwell’s equations using the
1820
FDFD method. In the vicinity of the waveguide junction,
higher order modes will be excited. We chose the left and
right ports of our junction sufficiently (5) away from the
physical junction where the amplitudes of the higher order
modes are negligible. In the following we will formulate
S in terms of the transverse magnetic field component Hy .
We will use &L x to denote the main mode of the left
waveguide and &R x for the right waveguide.
The scattering matrix relates the amplitudes and phases
of the modes that arrive at the left and right ports, HL+ , HR+
to the amplitudes and phases of the modes that propagate
away from the ports, HL− and HR− . Formally we can write
HL−
HR−
=
S
11
S21
S12 H +
L
S22
HR+
(37)
In order to extract S we need to know the fields that
arrive at the left and right ports from our numerical sources
in the simulation domain. To do that we do two calibration simulations (one for the left waveguide, another
for the right waveguide) without any discontinuities, as
shown in Figure 7(b), and record the fields. This gives us
+
the required H/L
R2 &/L R2 x in addition to the propagation vectors, k/L R2 , of the two main modes for guides of
insulator thicknesses 2a and 2a′ respectively. Then we do
two more simulations where we send the mode from the
left and from the right waveguide to the discontinuity as
schematically shown in Figures 7(c–e). From the results
of the simulation in Figure 7(c), for the fields to the left
of the left port, HL x z, we get
HL x z = HL+ e−ikL z + HL− e+ikL z &L x
= HL+ e−ikL z + S11 e+ikL z &L x
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Veronis et al.
Modeling of Plasmonic Waveguide Components and Networks
where the location of the left port determines the origin
for z and in (37) we used the fact that HR+ = 0 for the
simulation depicted in Figure 7(c). Simple algebra gives
Taking these circuit and transmission line approaches
together, we can then model a broad range of MIM systems with circuit models.
There is no unique way to describe S using lumped cirH x z −ikL z
cuit elements.71 To choose one circuit out of the infinite
− e−2ikL z
(38)
S11 = +L
e
HL &L x
possible set that could correspond to the same S, we will
first look at the well studied PEC case. The solution to
Very similarly, we also extract S22 from the results of the
the scattering problem for the junction of two PEC parsimulation of Figure 7(d).
allel plate waveguides was developed and experimentally
In order to extract S12 , we terminate our simulation
verified.72 It consists of a capacitor with susceptance B and
domain at the plane of the right port with a perfect electric
a transformer with a turns ratio of n:1. The susceptance
conductor. Such a termination results in zero tangential
and the turns ratio are described in terms of the geometry
electric fields and therefore gives −1 for the reflection
of the junction. The square of the turns ratio of the transcoefficient of the transverse electric field, Ex , and +1 for
former is equal to
the magnetic field, Hy . Thus, at the right port we get HR− =
a′
HR+ . Using this equality in (37), we get the reflection coefn2 =
(41)
a
ficient, R = HL− /HL+ , in Figure 7(e)
It is worthwhile remembering that the primary-secondary
2
S12
S12 S21
turns ratio of the transformer, n:1, is also the ratio of
= S11 +
(39)
R = S11 +
the voltages
1 − S22
1 − S22 Delivered by Ingenta
to: at its terminals. From the conservation of
power, currents have the inverse ratio and as a result
Houston
where in the last equality we used the fact thatUniversity
S21 = S12 . ofthe
impedance ratio at the transformer terminals is n2 :1.
We extract R using the same method as we usedIPin: 129.7.158.43
the
The derivation of the circuit elements for the PEC case
Mon,
05 Oct
14:12:14
S11 and
S22 2009
extraction of S11 . From the knowledge of R,
can be found in Refs. [73–75]. Note that for the PEC
one can easily invert (39) to calculate S12 .
case, only two parameters, B and n, are sufficient to
After we calculate S for the ports defined in Figure 7(a),
describe the junction even though in general three paramewe shift both the left and right reference planes back to
ters are required for a lossless reciprocal system. The nonthe exact location of the junction using
dispersive nature of the main mode of PEC parallel plate
ik ℓ
ik ℓ
waveguides leads to a further symmetry in the junction
L
L
L
L
0
e
0
e
(40)
SJ =
ikR ℓR
ikR ℓR S
which reduces the number of circuit parameters required.
0
e
0
e
At optical frequencies where the modes are strongly diswhere ℓL = ℓR = 5 as defined in Figure 7(a) and SJ is
persive, a third circuit element is needed in order to be
the effective scattering matrix for the case where the left
able to fit the elements of S exactly. For that reason we
and right ports are projected back to coincide with the
have an inductor term with a reactance X. A schematic of
junction plane.69 For the sake of notational abbreviation,
the circuit diagram is shown in Figure 8(a). The PEC parfrom this point on we will use S to imply SJ . Note that
allel plate waveguide circuit is the same, with X = 0. The
this effective scattering matrix is defined just for the algenormalization that we defined in (36) leads to transmission
braic convenience of having a scattering matrix associated
lines with a unit characteristic impedance on both sides
directly with the position of the interface. The fields near
of the junction. From transmission line theory we get the
to the interface are not in fact just describable by these
following equalities in terms of the equivalent impedance
single main modes because of various near field effects of
looking from the left side of the circuit, Z1 , and the equivhigher order modes.
alent admittance looking from the right side, Y2
Note that the polar plot of the reflection coefficients
Z −1
1 − Y2
is nothing other than the Smith chart of the microwave
−S11 = 1
and − S22 =
(42)
70
Z
+
1
1 + Y2
theory, which we will use in Section 5.4.
1
where
5.3. Circuit Model for the Waveguide Junction
Another important approach in microwave waveguide
modeling is the use of equivalent circuit models, which can
give an intuitive picture of the system as well as allowing
the use of circuit simulators for design.
Here we relate the scattering matrix to simplified circuit
models that can characterize the MIM waveguide interfaces. Since we only have single propagating modes in the
guides we consider, we can also use equivalent transmission lines to describe the propagation between interfaces.
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
1
iB + 1/n2
1
+ iB n2
Y2 =
1 + iX
Z1 = iX +
(43)
The reason why we have negative signs in front of S11 and
S22 in (42) is because we defined S based on the transverse
magnetic component of the main mode, Hy . However, the
norm in circuit parametrization is to use the voltage reflection and transmission coefficients, which correspond to
1821
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
1
0.8
iX
Y2
0.6
n2
iB
0.4
Z1
n:1
0.2
(a)
(b)
0
0
0.2
0.4
0.6
0.8
1
a′/a
0.07
2
(c)
(d)
0.06
1.6
0.05
1.2
X
B
0.04
Delivered by Ingenta
to:
University of Houston
0.03
IP : 129.7.158.43
Mon, 05 Oct 20090.0214:12:14
0.8
0.4
0.01
0
0
0
0.2
0.4
0.6
0.8
1
a′/a
0
0.2
0.4
0.6
0.8
1
a′/a
Fig. 8. (a) Circuit description of the lossless ( m = −143497) waveguide junction for 2a = 01 (), 2a = 05 (), 2a = 09 (•). (b) Square of
the turns ratio, n2 , which is equivalent to the impedance ratio at the terminals of the transformer. Dashed line is the PEC result described by (41).
(c) Susceptance, B, for the MIM (, , •) and the PEC (△, , ) case. (d) Reactance, X, of the MIM waveguide. X = 0 for the PEC parallel plate
waveguide.
a scattering matrix description for the transverse electric
component, Ex . Just as in transmission line theory where
the reflection coefficient for voltage is the negative of that
of the current, the same relationship also holds exactly
between the reflection coefficients of Ex and Hy .
Let the real and imaginary parts of Z1 and Y2 be denoted
as ZR = ReZ1 , ZI = ImZ1 , YR = ReY2 and YI =
ImY2 . Using (43) we get
Bn2 =
Y I + YR ZI
1 − YR ZR
and X =
Z I + Z R YI
1 − YR ZR
(44)
Once we know Bn2 and X, we can calculate n2 using
(43) as
n2 = YR 1 + X 2 = ZR 1 + Bn2 2
(45)
Using Eqs. (42), (44) and (45) one can calculate the circuit
parameters from S11 and S22 . In Figures 8(b–d) we plotted
n2 , B and X as a function of a′ /a for the three different
fixed 2a values of 01, 05 and 09. It can be seen that
the PEC circuit description and the MIM circuit description lead to parameters which qualitatively have similar
behaviors.
1822
5.4. Cascade Connection of Junctions
We now test the utility of the scattering matrix description
by numerically simulating mode propagation through a
cascade connection of junctions and comparing the results
with the predictions of the scattering matrix formalism.
When different scattering matrices are cascaded, the
transfer matrix, T, leads to a much simpler formulation76
+
−
HL
T11 T12
HR
(46)
=
HR+
T21 T22
HL−
5.4.1. Conditions for Zero Reflection
In order to have HL− = 0, one should have T21 HR− +
T22 HR+ = 0, which can be cast in terms of the scattering
parameters as
S11 HR− = S11 S22 − S12 S21 HR+
(47)
a. Lossy Case. Let us investigate the case when two
junctions characterized by two different scattering matrices, L S and R S, are separated by a center waveguide of
length ℓ as shown in Figure 9(a). Suppose that we adjust
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
–
+
HC
HL
(a)
RS
1
ℓ
LS
LS
LS
11
21
LS
RS
12
RS
22
–
+
HL
RS
HC
11
21
RS
RS
21
12
22
11
90º
(b)
(c)
0.7
PL
0.6
22
PR
180º
0º
RS
0.2
0.4
0.6
0.8
1.
270º
11
Power reflection
LS
0.5
0.4
0.3
0.2
Delivered by Ingenta to:
0.1
University of Houston
IP : 129.7.158.43
0
ℓ
0
0.2
300Mon, 05 Oct 2009 14:12:14
0.9λ 0.5λ
w
0.4
0.6
0.8
1
1.2
Center waveguide length (ℓ/λ)
Fig. 9. (a) Schematic diagram of modal propagation. The left and right junctions are shown as boxes with an S matrix description. The center
waveguide is shown as a transmission line of length ℓ. The source that creates the fields is normalized such that the mode that propagates to the right
has a unit magnitude at the input of the right junction. (b) Graphical solution of (50) and (51) on the complex plane. Point PL is the location of the
left junction on the L S22 curve where 2a = 09 and 2a′ = 05. Point PR is the location of the right junction on the R S11 curve where 2a = 05 and
2a′ = 016. (c) Test of the scattering matrix description. Horizontal axis is the length of the center waveguide normalized to . Vertical axis is the
power reflection coefficient. FDFD simulation results •, transfer matrix calculations using lossy junctions (solid line) and lossless junctions (dashed
line) are also plotted. Transfer matrix calculations do take into account the loss in the center waveguide for both cases. As the junctions get very close
to each other <01 transfer matrix model begins to break down due to higher order modal interactions.
our excitation amplitude such that the mode that propagates toward the right junction at its input plane, which is
the junction plane, has unit strength. That choice of normalization leads to HR− = eikC ℓ and HR+ = R S11 e−ikC ℓ where
kC is the wave vector of the center waveguide. With these
definitions, the condition for zero reflection, (47), for the
left junction can be written as
L
e−2ikC ℓ =
LS
11
S11 /R S11
L
L
22 − S12 S21
LS
For reciprocal media, S12 = S21 , we can write
L
e−2ikC ℓ =
S11 /R S11
LS LS − LS2
11
22
21
(48)
b. Lossless Case. If the system is lossless, then the scattering matrix should be unitary (SS† = 1) which implies
the following three conditions
S11 2 = S22 2 = 1 − S12 2
S
S12
= − 22
∗
∗
S21
S11
(49)
Using (48) and (49), the zero reflection condition becomes
L S22 = R S11
(magnitude condition)
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
(50)
∡L S22 + ∡R S11 = 2kC ℓ + 2;n (phase condition) (51)
where n is any integer value, superscripts R and L denote
right and left respectively. “∡” is used to represent the
argument of a complex number. What this means is that,
to match a left waveguide to a right waveguide, one should
choose a center waveguide width which satisfies the magnitude condition, and decide on the length of the center
waveguide based on the phase condition.
Once a matching left, center and right waveguide triplet
is found, the procedure can be recursively repeated to cascade more junctions without getting any reflection at the
leftmost waveguide.
5.4.2. Mode Converter Design
Now that we have the conditions (48) and (50–51) for
zero reflection, we can test their validity. Condition (48) is
more general and is applicable to the lossy case. We did a
series of simulations in which we extracted S for a hypothetical lossless metal with a real, negative permittivity
m = −143497. The results were very similar to the case
1823
Modeling of Plasmonic Waveguide Components and Networks
Veronis et al.
where the loss was included. That led us to suspect that
6. SUMMARY AND CONCLUSIONS
the conditions for the lossless reciprocal junctions, (50)
In this paper, we first gave an overview of the challenges
and (51), would be essentially sufficient in the design of
involved in modeling of plasmonic devices, and briefly
a mode converter that converts the mode of a wavelength
examined two general purpose simulation techniques which
sized MIM waveguide 2a = 09 to that of a subwaveare widely used for modeling of plasmonic waveguide
length one with no reflection.
devices.
In our design we choose the left waveguide width to be
We then introduced the mode-matching technique for
09 and the center waveguide width to be 05 as shown
modeling of MIM waveguide devices. We derived the
in the inset in Figure 9(b). The parameters that we need
complete set of orthogonal modes that the MIM wavegare the insulator width of the right waveguide, w, and the
uide supports. We showed that in addition to the previlength of the center waveguide, ℓ.
ously unreported complex set of discrete modes, there also
The width of the right waveguide can be chosen by satexists a continuous set of modes. We then used the comisfying (50). In Figure 9(b) PL is the location of the 2a =
plete set of modes to apply the mode-matching technique
09 to 2a′ = 05 junction on the L S22 curve. To satisfy
in the investigation of modal scattering at the symmet(50) we need to have L S22 = R S11 . The solution can be
ric junction of two MIM waveguides with different cross
graphically found by drawing a circle in the complex plane
sections.
with a radius PL centered at the origin and finding its
This special purpose simulation technique is far more
intersection with the R S11 curve. The intersection point is
efficient
Delivered by Ingenta for
to: this class of problems than general purpose
denoted by PR . PR corresponds to a right waveguide thickelectromagnetic
simulation techniques. The knowledge of
University of Houston
ness of 016. The phase condition (51) is then easily
calthe
set
of
orthogonal
modes which form a complete basis
IP : 129.7.158.43
culated from the phases of the scattering coefficients,
∡P
for
a
given
geometry
leads to a much more simplified
L
Mon, 05 Oct 2009 14:12:14
algebra and speeds up calculations. These results are valuand ∡PR . After some simplification through the use of the
able for electromagnetic scattering calculations involving
numerical value for kC one gets ℓ/ = 01377 + 04861n,
the MIM geometry. The analysis can also be generalized
where n is any positive integer.
for other related geometries involving metals at optical
To check our design, we numerically simulated the
frequencies.77–81
structure shown in the inset of Figure 9(b) using FDFD and
We also introduced several different equivalent modlooked at the amount of power reflected back as a function
els
for plasmonic waveguide components. We first showed
of the center waveguide length ℓ. We also calculated the
that
the characteristic impedance model can accurately
power reflection coefficient through the use of the transmodel
the behavior of deep subwavelength MIM wavegfer matrix formalism in which we multiplied the transfer
uides.
We
also extracted the scattering matrices of juncmatrices for the right junction, TR , a center waveguide of
tions
of
different
geometries from full field solutions, and
length ℓ, TC , and the left junction, TL , to get the overparametrized
the
scattering
matrix of the MIM junction in
2
all transfer matrix T = TL TC TR , and plotted T21 /T11 of
terms
of
lumped
circuit
elements.
We validated our charT as a function of ℓ/. We did the calculations for two
acterization by designing a mode converter that concendifferent sets of /TR TL 2: one in which we used the scattrates light from an MIM waveguide of wavelength-sized
tering matrices for the lossy junctions and another for the
dimension to one of subwavelength dimension with zero
lossless junctions. The center waveguide of length ℓ had
reflection.
loss in both cases i.e., kC = 103 − i945 × 10−4 2;/.
The model abstraction provided by these equivalent
Figure 9(c) verifies that lossless junction models are
models is important for the analysis and synthesis of
quite effective at modeling the waveguide discontinuities
device functions.82–85 The circuit representation of the
and the prediction of the length of the center guide for zero
junction helps us associate the effects of geometry, matereflection reached by their use, ℓ/ = 01377 + 04861n
rial properties and wave propagation in terms of a simis very accurate. The lossy junction model on the other
ple network of a capacitor, inductor and a transformer.
hand gives results essentially indistinguishable from the
The scattering matrix description of junctions can be
simulation results as long as the two junctions are not
used to design optical circuitry with complex functionvery close to each other (< 01). When the junctions get
ality using tools of circuit analysis.86 87 It is conceivvery close, the coupling of higher order non-propagating
able to build a library of junction geometries associated
modes becomes important and the single mode modeling
with their scattering matrices for different waveguides
we employed in the construction of scattering matrices
including three dimensional nano-metallic ones.88 Such a
breaks down. For such closely spaced junctions, the whole
library, indexed according to modal scattering and propstructure should be treated as a single unit and its characagation properties, would be invaluable in the design of
teristics should be extracted by the techniques described
integrated optical circuits composed of many interacting
components.
in Section 5.2.
1824
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
Veronis et al.
Modeling of Plasmonic Waveguide Components and Networks
References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40. I. Breukelaar, R. Charbonneau, and P. Berini, J. Appl. Phys. 100,
043104 (2006).
W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature 424, 824
41. R. F. Oulton, D. F. P. Pile, Y. Liu, and X. Zhang, Phys. Rev. B 76,
(2003).
035408 (2007).
E. Ozbay, Science 311, 189 (2006).
42. Ş. E. Kocabaş, G. Veronis, D. A. B. Miller, and S. Fan, Phys. Rev.
J. Takahara, S. Yamagishi, H. Taki, A. Morimoto, and T. Kobayashi,
B 79, 035120 (2009).
Opt. Lett. 22, 475 (1997).
43. G. W. Hanson and A. B. Yakovlev, Operator Theory for ElectromagJ. C. Weeber, A. Dereux, C. Girard, J. R. Krenn, and J. P. Goudonnet,
netics, Springer-Verlag, New York (2002).
Phys. Rev. B 60, 9061 (1999).
44. B. Prade, J. Y. Vinet, and A. Mysyrowicz, Phys. Rev. B 44, 13556
J. R. Krenn, B. Lamprecht, H. Ditlbacher, G. Schider, M. Salerno,
(1991).
A. Leitner, and F. R. Aussenegg, Europhys. Lett. 60, 663 (2002).
45. E. Feigenbaum and M. Orenstein, Journal of Lightwave Technology
M. L. Brongersma, J. W. Hartman, and H. A. Atwater, Phys. Rev. B
25, 2547 (2007b).
62, R16356 (2000).
46. A. R. Zakharian, J. V. Moloney, and M. Mansuripur, Opt. Express
S. A. Maier, P. G. Kik, H. A. Atwater, S. Meltzer, E. Harel, B. E.
15, 183 (2007).
Koel, and A. A. G. Requicha, Nat. Mater. 2, 229 (2003).
47. V. V. Shevchenko, Continuous Transitions in Open Waveguides,
S. I. Bozhevolnyi, V. S. Volkov, E. Devaux, and T. W. Ebbesen,
Golem Press, Boulder, Colorado (1971).
Phys. Rev. Lett. 95, 046802 (2005).
48.
T.
Rozzi and M. Mongiardo, Open Electromagnetic Waveguides, The
E. N. Economou, Phys. Rev. 182, 539 (1969).
Institution of Electrical Engineers, London (1997).
R. Zia, M. D. Selker, P. B. Catrysse, and M. L. Brongersma, J. Opt.
49. D. A. B. Miller, Quantum Mechanics for Scientists and Engineers,
Soc. Am. A 21, 2442 (2004).
Cambridge University Press, Cambridge (2008).
J. A. Dionne, L. A. Sweatlock, H. A. Atwater, and A. Polman, Phys.
50.
M. Mrozowski, Guided Electromagnetic Waves, Research Studies
Rev. B 73, 035407 (2006a).
Press, UK (1997).
E. Feigenbaum and M. Orenstein Journal of Lightwave
Technology
Delivered by Ingenta
to:
51. W. C. Chew, Waves and Fields in Inhomogenous Media, Van Nos25, 2547 (2007a).
University
of
Houston
trand Reinhold, New York (1990).
R. Gordon, Phys. Rev. B 73, 153405 (2006).
IP : 129.7.158.43
52. R. E. Collin, Field Theory of Guided Waves, 2nd edn., WileyZ. Han and S. He, Opt. Commun. 278, 199 (2007).
Mon,
05 Oct 2009 Interscience,
14:12:14New York (1991a).
G. Veronis and S. Fan, Appl. Phys. Lett. 87, 131102
(2005).
53. E. Anemogiannis and E. N. Glytsis, Journal of Lightwave TechnolK. Tanaka and M. Tanaka, Appl. Phys. Lett. 82, 1158 (2003).
ogy 10, 1344 (1992).
G. Veronis and S. Fan, Opt. Express 15, 1211 (2007a).
54. A. Bakhtazad, H. Abiri, and R. Ghayour, Journal of Lightwave TechP. Ginzburg and M. Orenstein, Opt. Express 15, 6762 (2007).
J. A. Dionne, H. J. Lezec, and H. A. Atwater, Nano Lett. 6, 1928
nology 15, 383 (1997).
(2006b).
55. M. S. Kwon and S. Y. Shin, Opt. Commun. 233, 119 (2004).
Y. Jiao, S. Fan, and D. Miller, IEEE Journal of Quantum Electronics
56. R. Rodriguez-Berral, F. Mesa, and F. Medina, IEEE Transactions on
42, 266 (2006).
Microwave Theory and Techniques 53, 1613 (2005).
W. Suh, Z. Wang, and S. Fan, IEEE Journal of Quantum Electronics
57. R. Rodriguez-Berral, F. Mesa, and F. Medina, International Journal
40, 1511 (2004).
of RF and Microwave Computer-Aided Engineering 14, 73 (2004).
J. Bravo-Abad, S. Fan, S. G. Johnson, J. D. Joannopoulos, and
58. R. E. Smith, S. N. Houde-Walter, and G. W. Forbes, IEEE Journal
M. Soljacic, Journal of Lightwave Technology 25, 2539 (2007).
of Quantum Electronics 28, 1520 (1992).
S. Fan, M. F. Yanik, Z. Wang, S. Sandhu, and M. L. Povinelli,
59. H. J. Hagemann, W. Gudat, and C. Kunz, Journal of the Optical
Journal of Lightwave Technology 24, 4493 (2006).
Society of America 65, 742 (1975).
L. Novotny, B. Hecht, and D. Pohl J. Appl. Phys. 81, 1798 (1997).
60. A. Zettl, Sturm-Liouville Theory, American Mathematical Society
E. Prodan, P. Nordlander, and N. Halas, Chem. Phys. Lett. 368, 94
(2005).
(2003).
61. L. P. Felsen and N. Marcuvitz, Radiation and Scattering of Waves,
C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light
IEEE Press, New York (1994).
by Small Particles, Wiley Science, New York (1983).
62. G. V. Eleftheriades, A. S. Omar, L. P. B. Katehi, and G. M. Rebeiz,
E. D. Palik, Handbook of Optical Constants of Solids, Academic,
IEEE Transactions on Microwave Theory and Techniques 42, 1896
New York (1985).
(1994).
A. Taflove and S. C. Hagness, Computational Electrodynamics,
63. A. A. Oliner, IEEE Transactions on Microwave Theory and TechArtech House, Boston (2005).
niques MTT/32, 1022 (1984).
A. D. Rakic, A. B. Djurisic, J. M. Elazar, and M. L. Majewski, Appl.
64. S. A. Schelkunoff, IEEE Transactions on Antennas and Propagation
Opt. 37, 5271 (1998).
AP/18, 309 (1970).
M. Han, R. Dutton, and S. Fan, IEEE Microwave and Wireless Com65. P. J. Schuck, D. P. Fromm, A. Sundaramurthy, G. S. Kino, and W. E.
ponents Letters 16, 119 (2006).
Moerner, Phys. Rev. Lett. 94, 017402 (2005).
A. Vial, A. Grimault, D. Macias, D. Barchiesi, and M. de la
66. T. Takano and J. Hamasaki, IEEE Journal of Quantum Electronics
Chapelle, Phys. Rev. B 71, 85416 (2005).
QE/8, 206 (1972).
J. Jin, The Finite Element Method in Electromagnetics, John Wiley
67. S. Ramo, J. R. Whinnery, and T. V. Duzer, Fields and Waves in
& Sons, New York (2002).
Communication Electronics, Wiley, New York (1994).
P. Berini, Phys. Rev. B 61, 10484 (2000).
68. Ş. E. Kocabaş, G. Veronis, D. A. B. Miller, and S. Fan, IEEE Journal
G. Veronis, R. W. Dutton, and S. Fan, Opt. Lett. 29, 2288 (2004).
of Selected Topics in Quantum Electronics 14, 1462 (2008b).
J. Berenger, J. Comput. Phys. 114, 185 (1994).
69. D. M. Pozar, Microwave Engineering, Wiley, New York (1998).
J. Young and R. Nelson, IEEE Antennas and Propagation Magazine
70. P. H. Smith, Electronic Applications of the Smith Chart in Waveg43, 61 (2001).
uide, Circuit and Component Analysis, Robert E. Krieger Publishing
P. J. B. Clarricoats and K. R. Slinn, Proceedings of the Institution
Company, Malabar, FL (1983).
of Electrical Engineers (1967), Vol. 114, p. 878.
71. E. L. Ginzton, Microwave Measurements, McGraw-Hill Book ComP. Bienstman, Ph.D. thesis, Ghent University, Belgium (2001).
pany Inc. (1957); Representation and Measurement of Microwave
I. Breukelaar and P. Berini, Journal of the Optical Society of
America A 23, 1971 (2006).
Circuits, pp. 313–345.
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009
1825
Modeling of Plasmonic Waveguide Components and Networks
72. N. Marcuvitz, Waveguide Handbook, Peter Peregrinus Ltd., London,
UK (1986), Chap. 5, pp. 307–309.
73. J. Schwinger and D. S. Saxon, Discontinuities in Waveguides,
Gordon and Breach Science Publishers Inc., New York (1968),
Chap. 5, pp. 99–124.
74. L. Lewin, Advanced Theory of Waveguides, Iliffe & Sons Ltd.
(1951), Chap. 5, pp. 98–106.
75. R. E. Collin, Field Theory of Guided Waves, 2nd edn., WileyInterscience (1991b), Chap. 8, pp. 581–588.
76. R. E. Collin, Foundations for Microwave Engineering, McGraw-Hill
Book Company (1966).
77. S. Collin, F. Pardo, and J. L. Pelouard, Opt. Express 15, 4310 (2007).
78. F. M. Kong, K. Li, B. I. Wu, H. Huang, H. S. Chen, and J. A. Kong,
Progress in Electromagnetics Research 76, 449 (2007a).
79. F. M. Kong, B. I. Wu, H. S. Chen, and J. A. Kong, Opt. Express 15,
12331 (2007b).
Veronis et al.
80. F. Lopez-Tejeira, S. G. Rodrigo, L. Martin-Moreno, F. J. GarciaVidal, E. Devaux, T. W. Ebbesen, J. R. Krenn, I. P. Radko, S. I.
Bozhevolnyi, M. U. Gonzalez, et al., Nat. Phys. 3, 324 (2007).
81. M. W. Vogel and D. K. Gramotnev, Phys. Lett. A 363, 507 (2007).
82. A. Alu and N. Engheta, Journal of the Optical Society of America
B 23, 571 (2006).
83. N. Engheta, Science 317, 1698 (2007).
84. A. Hosseini, H. Nejati, and Y. Massoud, Opt. Express 15, 15280
(2007).
85. A. Hosseini, H. Nejati, and Y. Massoud, Opt. Express 16, 1475
(2008).
86. SPICE http://en.wikipedia.org/wiki/SPICE.
87. HSPICE® Signal Integrity User Guide Synopsys, z-2007.03 edn.,
(2007), ch. S-parameter Modeling Using the S-element.
88. G. Veronis and S. Fan, Journal of Lightwave Technology 25, 2511
(2007b).
Received: 8 October 2008. Accepted: 22 October 2008.
Delivered by Ingenta to:
University of Houston
IP : 129.7.158.43
Mon, 05 Oct 2009 14:12:14
1826
J. Comput. Theor. Nanosci. 6, 1808–1826, 2009