We are IntechOpen, the world's leading publisher of Open Access books Built by scientists, for scientists

4,200+

Open access books available

116,000+

International authors and editors

125M+

Downloads

151 Countries delivered to Our authors are among the

Top 1%

most cited scientists

12.2% Contributors from top 500 universities

Selection of our books indexed in the Book Citation Index in Web of Science™ Core Collection (BKCI)

## Interested in publishing with us? Contact book.department@intechopen.com

Numbers displayed above are based on latest data collected. For more information visit www.intechopen.com

Contents

**Preface IX** 

Chapter 1 **Modelling at Nanoscale 3**  Paolo Di Sia

Xingyu Gao

Guy A. E. Vandenbosch

**Section 1 Modeling and Computational Methods for Plasmonics 1** 

**Super-Resolution Focusing and Nano-Waveguiding 49** 

Chapter 2 **Computational Electromagnetics in Plasmonics 23** 

Chapter 3 **Numerical Simulations of Surface Plasmons** 

**Section 2 Plasmonic Structures for Light Transmission, Focusing and Guiding 73** 

**Nanofocusing in Asymmetric Structures 99** 

Chapter 6 **Propagating Surface Plasmons and Dispersion Relations** 

**for Nanoscale Multilayer Metallic-Dielectric Films 135**  Henrique T. M. C. M. Baltar, Krystyna Drozdowicz-Tomsia

Chapter 4 **Surface Plasmons on Complex Symmetry Nanostructured Arrays 75**  Brian Ashall and Dominic Zerulla

> Valeria Lotito, Christian Hafner, Urs Sennhauser and Gian-Luca Bona

Chapter 7 **Light Transmission via Subwavelength** 

**Apertures in Metallic Thin Films 157**  V. A. G. Rivera, F. A. Ferri, O. B. Silva, F. W. A. Sobreira and E. Marega Jr.

Chapter 5 **Novel SNOM Probes Based on** 

and Ewa M. Goldys

## Contents

#### **Preface XI**


V. A. G. Rivera, F. A. Ferri, O. B. Silva, F. W. A. Sobreira and E. Marega Jr.


Contents VII

Chapter 19 **Plasmonic Rectenna for** 

Chapter 20 **Application of Surface Plasmon** 

**Efficient Conversion of Light into Electricity 469** 

Fuyi Chen, Jian Liu and Negash Alemu

**Polaritons in CMOS Digital Imaging 495**  Qin Chen, Xiaohua Shi, Yong Ma and Jin He

Odysseas Tsilipakos, Alexandros Pitilakis, Emmanouil E. Kriezis and Nikos Pleros

Chapter 21 **Merging Plasmonics and Silicon Photonics Towards Greener** 

**and High-Performance Computing Systems 523**  Sotirios Papaioannou, Konstantinos Vyrsokinos, Dimitrios Kalavrouziotis, Giannis Giannoulis, Dimitrios Apostolopoulos, Hercules Avramopoulos, Filimon Zacharatos, Karim Hassan, Jean-Claude Weeber, Laurent Markey, Alain Dereux, Ashwani Kumar, Sergey I. Bozhevolnyi, Alpaslan Suna, Oriol Gili de Villasante, Tolga Tekin, Michael Waldow,

**and Faster "Network-on-Chip" Solutions for Data Centers** 

	- **Section 4 Plasmonics Applications 401**

Chapter 19 **Plasmonic Rectenna for Efficient Conversion of Light into Electricity 469**  Fuyi Chen, Jian Liu and Negash Alemu

#### Chapter 20 **Application of Surface Plasmon Polaritons in CMOS Digital Imaging 495**  Qin Chen, Xiaohua Shi, Yong Ma and Jin He

VI Contents

Chapter 8 **Plasmonic Lenses 183**

Yongqi Fu, Jun Wang and Daohua Zhang

Joong Wook Lee and DaiSik Kim

**Section 3 Emerging Concepts with Plasmonics 251** 

Chapter 10 **Application of Surface Plasmon Resonance Based on a Metal Nanoparticle 253**  Amir Reza Sadrolhosseini, A. S. M. Noor

Chapter 11 **Localized Surface Plasmon Resonances: Noble Metal** 

V.A.G. Rivera, F.A. Ferri and E. Marega Jr.

Chapter 12 **Resonant Excitation of Plasmons in Bi-Gratings 313** 

Taikei Suyama, Akira Matsushima, Yoichi Okuno and Toyonori Matsuda

Chapter 13 **A Treatise on Magnetic Surface Polaritons:** 

Chapter 14 **Photoacoustic Based Surface Plasmon** 

Young Chul Jun

**Section 4 Plasmonics Applications 401** 

Chapter 16 **Plasmonic Conducting Polymers** 

Witold Jacak

**for Heavy Metal Sensing 403**

Majid Reayi and Afarin Bahrami

Chapter 17 **Innovative Exploitation of Grating-Coupled** 

G. Ruffato, G. Zacco and F. Romanato

Chapter 18 **Plasmon Mediated Energy Transport in PV Systems**

Mahnaz M. Abdi, Wan Mahmood Mat Yunus,

**Surface Plasmon Resonance for Sensing 419**

**with Photo-Active Surface Modified Metallically in Nano-Scale and in Metallic Nano-Chains 445** 

**and Experimental Realization 335** Yu-Hang Yang and Ta-Jen Yen

Chapter 15 **Electrically-Driven Active Plasmonic Devices 383**

**Nanoparticle Interaction with Rare-Earth Ions 283** 

**Theoretical Background, Numerical Verification** 

**Resonance Spectroscopy: An Investigation 359** K. Sathiyamoorthy, C. Vijayan and V.M. Murukeshan

and Mohd. Maarof Moksin

Chapter 9 **Fundamental Role of Periodicity and Geometric Shape to Resonant Terahertz Transmission 229** 

> Chapter 21 **Merging Plasmonics and Silicon Photonics Towards Greener and Faster "Network-on-Chip" Solutions for Data Centers and High-Performance Computing Systems 523**  Sotirios Papaioannou, Konstantinos Vyrsokinos, Dimitrios Kalavrouziotis, Giannis Giannoulis, Dimitrios Apostolopoulos, Hercules Avramopoulos, Filimon Zacharatos, Karim Hassan, Jean-Claude Weeber, Laurent Markey, Alain Dereux, Ashwani Kumar, Sergey I. Bozhevolnyi, Alpaslan Suna, Oriol Gili de Villasante, Tolga Tekin, Michael Waldow, Odysseas Tsilipakos, Alexandros Pitilakis, Emmanouil E. Kriezis and Nikos Pleros

Preface

chapters in this book.

Although plasmonic phenomena have been identified more than half century ago, recent accomplishments in this technological area have witnessed remarkable developments in both fundamental principles and quite various practical applications due to the current advancements in computational resources and nanoscale fabrication techniques, elevated needs for novel applications including photovoltaics, subwavelength imaging, optical interconnects, spectroscopy, microscopy, lithography, etc., and recent efforts of combination with conventionally and newly developed physics and techniques such as extraordinary light transmission, metamaterials,

The book attempts to encompass state-of-the-art research in modern plasmonic technologies, resulting in 21 excellent chapters, which has been sub-categorized into 4 sections; (1) Modeling and computational methods for plasmonics, (2) Plasmonic structures for light transmission, focusing, and guiding, (3) Emerging concepts with plasmonics, and (4) Plasmonics applications. Some chapters straddle the border of two

This book concerns itself with advanced recent research results on plasmonic technologies covering the aforementioned topics. Although it is a collection of specific technological issues, I strongly believe that the readers can obtain generous and overall ideas and knowledge of the state-of-the-art technologies in plasmonics. I hope the readers will be inspired to start or to improve further their own research and technologies and to expand potential applications, hopefully triggering new interdisciplinary research work. Lastly, special words of thanks should go to all the scientists and engineers who have devoted a great deal of time to writing excellent

**Ki Young Kim**

graphenes, nanoparticles, terahertz radiation, and many more.

or more sections due to the inherent nature of plasmonics research.

## Preface

Although plasmonic phenomena have been identified more than half century ago, recent accomplishments in this technological area have witnessed remarkable developments in both fundamental principles and quite various practical applications due to the current advancements in computational resources and nanoscale fabrication techniques, elevated needs for novel applications including photovoltaics, subwavelength imaging, optical interconnects, spectroscopy, microscopy, lithography, etc., and recent efforts of combination with conventionally and newly developed physics and techniques such as extraordinary light transmission, metamaterials, graphenes, nanoparticles, terahertz radiation, and many more.

The book attempts to encompass state-of-the-art research in modern plasmonic technologies, resulting in 21 excellent chapters, which has been sub-categorized into 4 sections; (1) Modeling and computational methods for plasmonics, (2) Plasmonic structures for light transmission, focusing, and guiding, (3) Emerging concepts with plasmonics, and (4) Plasmonics applications. Some chapters straddle the border of two or more sections due to the inherent nature of plasmonics research.

This book concerns itself with advanced recent research results on plasmonic technologies covering the aforementioned topics. Although it is a collection of specific technological issues, I strongly believe that the readers can obtain generous and overall ideas and knowledge of the state-of-the-art technologies in plasmonics. I hope the readers will be inspired to start or to improve further their own research and technologies and to expand potential applications, hopefully triggering new interdisciplinary research work. Lastly, special words of thanks should go to all the scientists and engineers who have devoted a great deal of time to writing excellent chapters in this book.

**Ki Young Kim**

**Section 1** 

**Modeling and Computational Methods** 

**for Plasmonics** 

**Modeling and Computational Methods for Plasmonics** 

**Chapter 1** 

© 2012 Di Sia, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

© 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

distribution, and reproduction in any medium, provided the original work is properly cited.

and reproduction in any medium, provided the original work is properly cited.

One of the most important aspects at nanoscale concerns the charge transport, which can be influenced by particles dimensions and assumes different characteristics with respect to those of bulk. In particular, if the mean free path of charges, due to scattering phenomena, is larger than the particle dimensions, we have a mesoscopic system, in which the transport depends on dimensions and in principle it is possible to correct the transport bulk theories by considering this phenomenon. A similar situation occurs also in a thin film, in which the smallest nanostructure dimension can be less than the free displacement and therefore requires variations to existing theoretical transport bulk models. Therefore a rigorous knowledge of transport properties is to be acquired. From a theoretical viewpoint, various techniques can be used for the comprehension of transport phenomena, in particular analytical descriptions based on transport equations. Many existing and used theories at today regard numerical approaches, not offering analytical results, which would be of great mathematical interest and in every case suitable to be implemented through the experimental data existing in literature and continuously found by the experimentalists.

To establish the applicability limit of a bulk model and to investigate the time response of systems at nanoscale, recently it has appeared a novel theoretical approach, based on correlation functions obtained by a complete Fourier transform of the frequency-dependent conductivity of the studied system [1]. With this model it is possible to calculate exactly the expressions of the most important connected functions, i.e. the velocities correlation

At nanoscale we are in the middle between classical and quantum effects, macroscopic and microscopic properties of the matter [2]. Actually one of the most important experimental technique for the study of the frequency-dependent complex-valued far-infrared photoconductivity σ ω( ) is the Time-resolved THz Spectroscopy (TRTS), an ultrafast noncontact optical probe; data are normally fitted via Drude-Lorentz, Drude-Smith or effective medium models [3]. The quantum effects in the nanoworld have opened new ways in a lot

function, the mean free displacement and the diffusion coefficient.

of old and new technological sectors.

**Modelling at Nanoscale** 

Additional information is available at the end of the chapter

Paolo Di Sia

http://dx.doi.org/10.5772/50755

**1. Introduction** 

### **Chapter 1**

## **Modelling at Nanoscale**

Paolo Di Sia

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50755

### **1. Introduction**

One of the most important aspects at nanoscale concerns the charge transport, which can be influenced by particles dimensions and assumes different characteristics with respect to those of bulk. In particular, if the mean free path of charges, due to scattering phenomena, is larger than the particle dimensions, we have a mesoscopic system, in which the transport depends on dimensions and in principle it is possible to correct the transport bulk theories by considering this phenomenon. A similar situation occurs also in a thin film, in which the smallest nanostructure dimension can be less than the free displacement and therefore requires variations to existing theoretical transport bulk models. Therefore a rigorous knowledge of transport properties is to be acquired. From a theoretical viewpoint, various techniques can be used for the comprehension of transport phenomena, in particular analytical descriptions based on transport equations. Many existing and used theories at today regard numerical approaches, not offering analytical results, which would be of great mathematical interest and in every case suitable to be implemented through the experimental data existing in literature and continuously found by the experimentalists.

To establish the applicability limit of a bulk model and to investigate the time response of systems at nanoscale, recently it has appeared a novel theoretical approach, based on correlation functions obtained by a complete Fourier transform of the frequency-dependent conductivity of the studied system [1]. With this model it is possible to calculate exactly the expressions of the most important connected functions, i.e. the velocities correlation function, the mean free displacement and the diffusion coefficient.

At nanoscale we are in the middle between classical and quantum effects, macroscopic and microscopic properties of the matter [2]. Actually one of the most important experimental technique for the study of the frequency-dependent complex-valued far-infrared photoconductivity σ ω( ) is the Time-resolved THz Spectroscopy (TRTS), an ultrafast noncontact optical probe; data are normally fitted via Drude-Lorentz, Drude-Smith or effective medium models [3]. The quantum effects in the nanoworld have opened new ways in a lot of old and new technological sectors.

© 2012 Di Sia, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 4 Plasmonics – Principles and Applications

In the following it is illustrated an overview of the fundamental models that, starting from the Drude model, have attempted to analyze and explain in increasing accuracy the transport phenomena of the matter at solid state level and in particular at nanoscale, until the recent developments concerning the variations of Drude-Lorentz-like models, which involve in particular the concept of plasmon.

Modelling at Nanoscale 5

(1)

, with

proportional to the applied field 0 *J E* =σ

corresponds to a current density *J*

has also underlined series difficulties.

**4. The most utilized Drude-Lorentz-like models** 

leading to a dielectric function of the form:

dispersed in a surrounding host matrix.

(BR) model", which is a variation of the MG model.

*m*

examine directly the possible compatibility with the Einstein relation.

2

following:

[8]:

phonons.

<sup>0</sup> σ= τ *ne m*/ (*n* is the electron density). Estimating the relaxation time *τ*, Drude and Lorentz have obtained values of conductivity in good accordance with the experiments. In presence of an electric field of the form 0 ( ) *i t Et E e*− ω <sup>=</sup> , the complex conductivity assumes the form <sup>0</sup> /(1 ) *i* <sup>ω</sup> σ =σ − ωτ . Such model, said "Drude-Lorentz model", has received great success, but

One of the most utilized models for describing experimental transport data is the Drude-Lorentz model [4,5]; with such model the main transport parameters are obtainable. Starting from the Drude-Lorentz model, it is possible to obtain the velocities correlation function, from this the quadratic average distance crossed by the charges as a function of time and

Considerable variations of this model were made in this years; the most used are the

1. the "Maxwell-Garnett (MG) model": in this model the dielectric function is given by a Drude term with an additional "vibrational" contribution at a finite frequency ω*o* [6],

// 2 2

()1 ( /)

ε ω= − +

2 2

*p s i i* ω ω

ω ω+ τ ω −ω − γω

where the amplitude ω*<sup>s</sup>* , the resonant frequency ω0 and the damping constant γ are material dependent constants. The MG model usually describes an isotropic matrix containing spherical inclusions isolated from each other, such as the metal particles

2. In the "effective medium theories (EMTs)" the electromagnetic interactions between pure materials and host matrices are approximately taken into account [7]. The commonly used EMTs include the Maxwell-Garnett (MG) model and the **"**Bruggeman

Generally, in the THz regime, the dielectric function ( ) *<sup>m</sup>*ε ω consists of contributions of the high-frequency dielectric constant, conduction free electrons, and lattice vibration

<sup>2</sup> 2 2 ( ) *j j*

ε ω =ε − +

2 2

*p st TO*

*<sup>j</sup> TO <sup>j</sup> i i* <sup>∞</sup> ω ε ω

in which ∞ε is the high-frequency dielectric constant, the second term describes the contribution of free electrons or plasmons and the last term stands for the optical

*j*

ω + γω ω −ω − Γ ω (2)

0

### **2. The Drude model**

The most meaningful characteristics of metals are their elevated properties of electric and thermal conduction. Over the years this fact has brought to think in terms of a model in which the electrons are relatively free and can move under the influence of electric fields. Historically two models of the elementary theory of metals are born:

The Drude model, published in 1900 and based on the kinetic theory of an electron gas in a solid. It is assumed that all the electrons have the same average kinetic energy *Em*;

The previous model integrated with the foundations of quantum mechanics, called Sommerfeld model.

In the Drude model the valence electrons of atoms are considered free inside the metal; all the electrons move as an electronic gas. It is assumed that all the electrons have the same energy *Em* and that it exists a mechanism of collisions among ions and electrons, allowing the thermal equilibrium for the electrons; this fact implies the application of the kinetic theory of gases to such electronic gas. Drude published the theory three years after the discovery of the electron from J. J. Thomson. The free electrons have only kinetic energy, not potential energy; the average energy *Em* is therefore (3/2)*kBT*. It can be correlated to an average quadratic speed *vm* from the relation *Em* = (3/2)*kBT* = *mvm*2/2, denoting *m* the mass of the free electron. At environment temperature it is *vm* ≅ 107 cm/s, representing the average thermal speed of the electrons. It was assumed also that the electrons have collisions as instant events, i.e. the time of the diffraction is very smaller with respect to every other considered time. Through such collisions, the electrons acquire a thermal equilibrium corresponding to the temperature *T* of the metal. If the electrons don´t collide, they move in linear way following the Newton laws. The presence of a constant electric field determines an extra average velocity (the drift velocity) given by *vd* = - (*eE*/*m*)*t*. The relaxation time *τ* is defined as the average time between two collisions; it is so possible to get a mean free path, defined by *lmfp* = *vm <sup>τ</sup>*. The current density is *cond <sup>J</sup>* =σ *<sup>E</sup>* , with σ*cond* electric conductivity. This result, obtained by Drude, is an important goal of the classic theory for the conduction of metals (said "Drude theory").

### **3. The Drude-Lorentz model**

The Lorentz model (1905) is a refining of the Drude model, in which the statistical aspects are specified. The electrons are considered as free charges, with charge "-*e*"; they are described by a maxwellian velocity distribution. Considering an electron gas in a spatial region with a constant electric field, the drift velocity of the electrons is constant; this corresponds to a current density *J* proportional to the applied field 0 *J E* =σ , with 2 <sup>0</sup> σ= τ *ne m*/ (*n* is the electron density). Estimating the relaxation time *τ*, Drude and Lorentz have obtained values of conductivity in good accordance with the experiments. In presence of an electric field of the form 0 ( ) *i t Et E e*− ω <sup>=</sup> , the complex conductivity assumes the form <sup>0</sup> /(1 ) *i* <sup>ω</sup> σ =σ − ωτ . Such model, said "Drude-Lorentz model", has received great success, but has also underlined series difficulties.

### **4. The most utilized Drude-Lorentz-like models**

4 Plasmonics – Principles and Applications

**2. The Drude model** 

Sommerfeld model.

involve in particular the concept of plasmon.

In the following it is illustrated an overview of the fundamental models that, starting from the Drude model, have attempted to analyze and explain in increasing accuracy the transport phenomena of the matter at solid state level and in particular at nanoscale, until the recent developments concerning the variations of Drude-Lorentz-like models, which

The most meaningful characteristics of metals are their elevated properties of electric and thermal conduction. Over the years this fact has brought to think in terms of a model in which the electrons are relatively free and can move under the influence of electric fields.

The Drude model, published in 1900 and based on the kinetic theory of an electron gas in a

The previous model integrated with the foundations of quantum mechanics, called

In the Drude model the valence electrons of atoms are considered free inside the metal; all the electrons move as an electronic gas. It is assumed that all the electrons have the same energy *Em* and that it exists a mechanism of collisions among ions and electrons, allowing the thermal equilibrium for the electrons; this fact implies the application of the kinetic theory of gases to such electronic gas. Drude published the theory three years after the discovery of the electron from J. J. Thomson. The free electrons have only kinetic energy, not potential energy; the average energy *Em* is therefore (3/2)*kBT*. It can be correlated to an average quadratic speed *vm* from the relation *Em* = (3/2)*kBT* = *mvm*2/2, denoting *m* the mass of the free electron. At environment temperature it is *vm* ≅ 107 cm/s, representing the average thermal speed of the electrons. It was assumed also that the electrons have collisions as instant events, i.e. the time of the diffraction is very smaller with respect to every other considered time. Through such collisions, the electrons acquire a thermal equilibrium corresponding to the temperature *T* of the metal. If the electrons don´t collide, they move in linear way following the Newton laws. The presence of a constant electric field determines an extra average velocity (the drift velocity) given by *vd* = - (*eE*/*m*)*t*. The relaxation time *τ* is defined as the average time between two collisions; it is so possible to get a mean free path,

result, obtained by Drude, is an important goal of the classic theory for the conduction of

The Lorentz model (1905) is a refining of the Drude model, in which the statistical aspects are specified. The electrons are considered as free charges, with charge "-*e*"; they are described by a maxwellian velocity distribution. Considering an electron gas in a spatial region with a constant electric field, the drift velocity of the electrons is constant; this

, with σ

*cond* electric conductivity. This

solid. It is assumed that all the electrons have the same average kinetic energy *Em*;

Historically two models of the elementary theory of metals are born:

defined by *lmfp* = *vm <sup>τ</sup>*. The current density is *cond <sup>J</sup>* =σ *<sup>E</sup>*

metals (said "Drude theory").

**3. The Drude-Lorentz model** 

One of the most utilized models for describing experimental transport data is the Drude-Lorentz model [4,5]; with such model the main transport parameters are obtainable. Starting from the Drude-Lorentz model, it is possible to obtain the velocities correlation function, from this the quadratic average distance crossed by the charges as a function of time and examine directly the possible compatibility with the Einstein relation.

Considerable variations of this model were made in this years; the most used are the following:

1. the "Maxwell-Garnett (MG) model": in this model the dielectric function is given by a Drude term with an additional "vibrational" contribution at a finite frequency ω*o* [6], leading to a dielectric function of the form:

$$\text{cc}\_{//}(\text{co}) = 1 - \frac{\text{co}\_p^2}{\text{co}(\text{co} + \text{i} / \text{ }\text{\textdegree\text})} + \frac{\text{co}\_s^2}{\text{co}\_0^2 - \text{co}^2 - \text{i} \,\text{\textdegree\textalpha}} \tag{1}$$

where the amplitude ω*<sup>s</sup>* , the resonant frequency ω0 and the damping constant γ are material dependent constants. The MG model usually describes an isotropic matrix containing spherical inclusions isolated from each other, such as the metal particles dispersed in a surrounding host matrix.

2. In the "effective medium theories (EMTs)" the electromagnetic interactions between pure materials and host matrices are approximately taken into account [7]. The commonly used EMTs include the Maxwell-Garnett (MG) model and the **"**Bruggeman (BR) model", which is a variation of the MG model.

Generally, in the THz regime, the dielectric function ( ) *<sup>m</sup>*ε ω consists of contributions of the high-frequency dielectric constant, conduction free electrons, and lattice vibration [8]:

$$\text{i.e.}\_{m}(\text{co}) = \text{e}\_{\text{co}} - \frac{\text{co}\_{p}^{2}}{\text{co}^{2} + i \text{γco}} + \sum\_{j} \frac{\text{e}\_{st\_{j}} \left| \text{o}\_{TO\_{j}} \right|^{2}}{\text{o}\_{TO\_{j}} \text{a}^{2} - \text{o}^{2} - i \text{ } \Gamma\_{j} \text{ } \text{co}} \tag{2}$$

in which ∞ε is the high-frequency dielectric constant, the second term describes the contribution of free electrons or plasmons and the last term stands for the optical phonons.

When the response originates mainly from the contribution of free electrons or plasmons, it is usually adopted the Drude model:

$$\mathfrak{e}\_m(\alpha) = \mathfrak{e}\_{\alpha} - \frac{\alpha \mathfrak{o}\_p^2}{\alpha^2 + i \gamma \mathfrak{o}\_p} \tag{3}$$

Modelling at Nanoscale 7

( )

*n <sup>n</sup> <sup>n</sup>*

(10)

(11)

1

=

1

τ

=

− ωτ − ωτ

∞

(1 ) (1 )

*n e c m i i*

σω= +

A lot of nanostructures, in particular the nanostructured metals, show very complex and interesting optical properties. One of the most interesting and promising phenomena encountered in these structures are electromagnetic resonances, due to collective oscillations of the conduction electrons, said "plasmons". Plasmon modes exist in different geometries and in various metals, in particular in noble metals such as gold and silver. Under determined circumstances, plasmons are excited by light; this leads to strong light scattering and absorption and to an increase of the local electromagnetic field. The interest in plasmon modes has started to the beginning of the 20th century, with Zenneck (1907), Mie (1908), Sommerfeld (1909) and other scientists; recent advances regarding the structure, the manipulation and the observation at nanoscale, have increased the study of this scientific topic. The theoretical efforts have encountered also the growing demand at technological level, in particular for semiconductor based integrated electronic components, optical applications and new nano-components. Actually it remains an important challenge for research and development processes, like the guide of light in integrated optical systems and the interface with electronic components. A lot of nanostructures are believed to be one

The concept of plasma is very useful in the description of some aspects of the interaction radiation/conductive matter; the free electrons of the conductive material (for example a metal) are considered as an electron fluid with high density, of order of <sup>23</sup> <sup>3</sup> 10 *cm*<sup>−</sup> . Such concept is the base of the classic Drude model, which assumes that the material contains stopped positive ions and a gas of classical not interacting electrons, whose motion results slowed by a force of viscous friction, due to the collisions with the ions, and characterized from a relaxation time τ . The motion of electrons results so casual, due to the continuous collisions with the lattice. The plasma frequency is defined as the proper frequency of the

∞

= −τ + τ

*n n*

( ) (0) exp( / ) 1 / ! *<sup>n</sup>*

The *nc* factors hold into account of the original electrons speed remained after the n-th

*jt j t c t n*

Eq. (9) can be considered the first term of a series of the form:

of the key ingredients of the future optoelectronic devices [10].

collective motion of electrons in the following way:

**6. Plasmonics** 

**6.1. Introduction** 

**6.2. Plasmons** 

collision. The analytical form of the complex conductivity is therefore [9]:

2

\* ( ) <sup>1</sup>

which described well the dielectric properties of metals and semiconductors [7].

If instead the interaction of a radiation field with the fundamental lattice vibration plays a dominant role and results in absorption of electromagnetic wave, due to the creation or annihilation of lattice vibration, the dielectric function ( ) *<sup>m</sup>*ε ω mainly consists of the contributions of the lattice vibrations, expressed by the classical pseudo-harmonic phonon model in the first approximation [7]:

$$\mathfrak{e}\_{m}(\mathsf{oo}) = \mathfrak{e}\_{\mathsf{oo}} + \frac{\mathfrak{e}\_{\mathrm{st}} \left\| \mathfrak{a}\_{\mathrm{TO}} \right\|^{2}}{\left\| \mathfrak{a}\_{\mathrm{TO}} \right\|^{2} - \mathfrak{a}^{2} - i \,\gamma \,\mathrm{co}} \tag{4}$$

#### **5. The Smith model**

Smith has started from the response theory for the optical conductivity, considering an electric field impulse applied to a system, in order to examine the answer with respect to the current. He has utilized the following Fourier transform for the frequency-dependent complex conductivity:

$$\sigma(\text{oi}) = \int\_0^\infty j(t) \exp(i\alpha t) dt \tag{5}$$

A field impulse, which exceeds every other acting force on the system, permits to assume the hypothesis that the electrons initially can be considered totally free; therefore it holds:

$$j(0) = n \, ^\ast e^2 \Big/ m \tag{6}$$

with *n*\* effective electron density. After calculation, the real part of σ ω( ) results:

$$\int\_0^\infty \text{Re}\,\sigma(\text{co})d\text{co} = \frac{\pi}{2}j(0) = \frac{\text{co}\_p^2}{8} \tag{7}$$

If the initial current decays exponentially to its initial value with relaxation time τ , it is possible to write:

$$(j(t)/j(0) = \exp(-t/\tau) \tag{8}$$

from which the standard Drude formula is obtainable:

$$\sigma(\alpha) = (n \, ^\ast e^2 \, ^\ast \pi / m) / (1 - i \, \alpha \, \tau) \tag{9}$$

Eq. (9) can be considered the first term of a series of the form:

$$j(t)\left/j(0) = \exp(-t/\tau)\left[1 + \sum\_{n=1}^{\infty} c\_n \left(t/\tau\right)^n / n!\right] \tag{10}$$

The *nc* factors hold into account of the original electrons speed remained after the n-th collision. The analytical form of the complex conductivity is therefore [9]:

$$\sigma(\alpha) = \frac{n^\* e^2 \pi}{m \left(1 - i \cot \mathfrak{r} \right)} \left[ 1 + \sum\_{n=1}^{\infty} \frac{c\_n}{\left(1 - i \cot \mathfrak{r} \right)^n} \right] \tag{11}$$

### **6. Plasmonics**

6 Plasmonics – Principles and Applications

is usually adopted the Drude model:

model in the first approximation [7]:

**5. The Smith model** 

complex conductivity:

possible to write:

When the response originates mainly from the contribution of free electrons or plasmons, it

<sup>2</sup> ( ) *<sup>p</sup> <sup>m</sup> i* <sup>∞</sup>

If instead the interaction of a radiation field with the fundamental lattice vibration plays a dominant role and results in absorption of electromagnetic wave, due to the creation or annihilation of lattice vibration, the dielectric function ( ) *<sup>m</sup>*ε ω mainly consists of the contributions of the lattice vibrations, expressed by the classical pseudo-harmonic phonon

2 2 ( ) *st TO*

Smith has started from the response theory for the optical conductivity, considering an electric field impulse applied to a system, in order to examine the answer with respect to the current. He has utilized the following Fourier transform for the frequency-dependent

<sup>0</sup> ( ) ( )exp( ) *j t i t dt*

A field impulse, which exceeds every other acting force on the system, permits to assume the hypothesis that the electrons initially can be considered totally free; therefore it holds:

> Re ( ) (0) 2 8 *<sup>p</sup> d j*

If the initial current decays exponentially to its initial value with relaxation time τ , it is

<sup>∞</sup> π ω

∞

with *n*\* effective electron density. After calculation, the real part of σ ω( ) results:

0

from which the standard Drude formula is obtainable:

*TO <sup>i</sup>* <sup>∞</sup> ε ω

*m*

ε ω =ε +

ε ω =ε −

which described well the dielectric properties of metals and semiconductors [7].

2

2

ω + γω (3)

ω −ω − γω (4)

σω= ω (5)

<sup>2</sup> *j ne m* (0) \* = (6)

σ ω ω= = (7)

*jt j t* ( ) (0) exp( / ) = −τ (8)

<sup>2</sup> σ ω = τ − ωτ ( ) ( \* / )/(1 ) *ne m i* (9)

2

ω

#### **6.1. Introduction**

A lot of nanostructures, in particular the nanostructured metals, show very complex and interesting optical properties. One of the most interesting and promising phenomena encountered in these structures are electromagnetic resonances, due to collective oscillations of the conduction electrons, said "plasmons". Plasmon modes exist in different geometries and in various metals, in particular in noble metals such as gold and silver. Under determined circumstances, plasmons are excited by light; this leads to strong light scattering and absorption and to an increase of the local electromagnetic field. The interest in plasmon modes has started to the beginning of the 20th century, with Zenneck (1907), Mie (1908), Sommerfeld (1909) and other scientists; recent advances regarding the structure, the manipulation and the observation at nanoscale, have increased the study of this scientific topic. The theoretical efforts have encountered also the growing demand at technological level, in particular for semiconductor based integrated electronic components, optical applications and new nano-components. Actually it remains an important challenge for research and development processes, like the guide of light in integrated optical systems and the interface with electronic components. A lot of nanostructures are believed to be one of the key ingredients of the future optoelectronic devices [10].

#### **6.2. Plasmons**

The concept of plasma is very useful in the description of some aspects of the interaction radiation/conductive matter; the free electrons of the conductive material (for example a metal) are considered as an electron fluid with high density, of order of <sup>23</sup> <sup>3</sup> 10 *cm*<sup>−</sup> . Such concept is the base of the classic Drude model, which assumes that the material contains stopped positive ions and a gas of classical not interacting electrons, whose motion results slowed by a force of viscous friction, due to the collisions with the ions, and characterized from a relaxation time τ . The motion of electrons results so casual, due to the continuous collisions with the lattice. The plasma frequency is defined as the proper frequency of the collective motion of electrons in the following way:

#### 8 Plasmonics – Principles and Applications

$$
\cos^2 \frac{\mathbf{o}\_0}{\mathbf{e}\_0 \cdot \mathbf{r}} = \frac{ne^2}{\mathbf{e}\_0 \cdot m} \tag{12}
$$

Modelling at Nanoscale 9

(14)

incident light, such theory predicts that the extinction caused by a metallic nanosphere is estimated in the quasi-static and dipole limit. The relative cross section of the process

3 3 2

ln(10) ( ) *Am i*

σ = <sup>λ</sup> ε +χε +ε

where *NA* is the surface density, *a* the radius of the metallic nanosphere, *m*ε the dielectric constant of the medium surrounding the nanosphere, λ the wavelength, *<sup>i</sup>* ε and *<sup>r</sup>* ε imaginary and real part of the dielectric function of the metallic nanosphere respectively and χ form factor concerning the nanoparticle. The localized resonance depends also on the

The Gans theory extends the Mie theory to spheroidal particles case. In the case of polarization of incident light parallel to the symmetry axis of the spheroid, the cross section

3 2 2

*m j j*

π ε ε

*j j*

( )

*P P*

( ) \*

*r mj j*

2 (1 / ) 3 1

*V P*

σ = <sup>λ</sup> <sup>−</sup> ε + ε +ε

with *V* volume of nanoparticle and *Pj* depolarization factors along the three cartesian aces,

With this method the nanoparticles are divided in a cubic array of *N* polarizable dipoles, with polarizabilities *<sup>i</sup>* α determined by the dielectric function of nanoparticles. The induced dipole *<sup>i</sup> P* of every element results *i i loc i*, *P E* =α in an applied plane wave field. The cross

> 2 , 1 <sup>4</sup> Im *N ext inc j j*

with *k c* =ω (in vacuum). Such method can be applied for the calculation of the absorption, scattering, extinction and other optical properties of nanoparticles of various forms and dimensions. It is considered one of the most important methods for the understanding of the structural characteristics and optical properties of nanomaterials

*<sup>k</sup> E P*

*<sup>j</sup> inc*

*E* <sup>=</sup> π σ = <sup>⋅</sup>

*N a* πε ε

24

*ext*

interparticle space and on the dielectric constant of the substrate.

*ext*

which hold into account of the anisotropic form of particles.

*6.3.3. The discrete-dipole approximation method (DDA)* 

2 2

2 2

(15)

(16)

*rm i*

results:

*6.3.2. The Gans theory* 

of the process is given by:

section is determined by:

and nanostructures.

with 0 σ conductivity of material and common meaning of the other mathematical symbols. The dielectric constant of metal can be written as a function to the plasma frequency; in first approximation it results:

$$
\varepsilon(\alpha) = 1 - \frac{\alpha\_p^2}{\alpha^2} \tag{13}
$$

From the comparison between ω and ω*<sup>p</sup>* , we can deduce the behaviour of the electromagnetic waves arriving to the metal. If it holds ω<ω*<sup>p</sup>* , the dielectric constant is negative, therefore its square root is imaginary pure; this involves the reflection of the incident wave. In the contrary case, the square root of the dielectric constant is real and the incident wave can propagate in the medium with a small attenuation.

The plasma oscillations are the fluctuations of charge density, which happen to the frequency ω*p* and propagate in the metal. The quanta of such fluctuations inside the volume are the "volume plasmons".

There is also a mechanism on the metal surface, characterized by the "surface plasmons". Normally it appears in the wavelength range between the visible and the infrared for the interface air/metal. Localized surface electromagnetic excitations can exist at the surface of nanoparticles and metallic nanostructures: the "localized surface plasmons" (LSP) [11]. The frequency and the intensity of the radiation are very sensitive with respect to the dimension, form and morphology of the nanostructures [12]. Nanoparticles and nanostructures with smaller dimensions with respect to the wavelength of the exciting light are characterized by a wide absorption band, normally in the range of the visible and near infrared spectrum.

### **6.3. Related theoretical models**

The models concerning plasmonics simulate the extinction of the "localized surface plasmons resonance" (LSPR) from nanoparticles and nanoarrays; the application regards the calculation of the light absorption and the scattering. The most used models are:


#### *6.3.1. The Mie theory*

The Mie theory is a theoretical approach concerning the optical properties of the nanoparticles. When the nanoparticle dimension is smaller than the wavelength of the incident light, such theory predicts that the extinction caused by a metallic nanosphere is estimated in the quasi-static and dipole limit. The relative cross section of the process results:

$$\sigma\_{ext} = \frac{24\pi N\_A a^3 \varepsilon\_m^{3/2}}{\lambda \ln(10)} \left[ \frac{\varepsilon\_i}{\left(\varepsilon\_r + \chi \varepsilon\_m\right)^2 + \varepsilon\_i^2} \right] \tag{14}$$

where *NA* is the surface density, *a* the radius of the metallic nanosphere, *m*ε the dielectric constant of the medium surrounding the nanosphere, λ the wavelength, *<sup>i</sup>* ε and *<sup>r</sup>* ε imaginary and real part of the dielectric function of the metallic nanosphere respectively and χ form factor concerning the nanoparticle. The localized resonance depends also on the interparticle space and on the dielectric constant of the substrate.

#### *6.3.2. The Gans theory*

8 Plasmonics – Principles and Applications

approximation it results:

are the "volume plasmons".

near infrared spectrum.

• the Mie theory [13]; • the Gans theory [14];

*6.3.1. The Mie theory* 

**6.3. Related theoretical models** 

2

(12)

(13)

*ne m*

2

ω

2 0

σ ω= = ετ ε

*p*

0 0

with 0 σ conductivity of material and common meaning of the other mathematical symbols. The dielectric constant of metal can be written as a function to the plasma frequency; in first

> <sup>2</sup> ()1 <sup>ω</sup>*<sup>p</sup>* εω= −

From the comparison between ω and ω*<sup>p</sup>* , we can deduce the behaviour of the electromagnetic waves arriving to the metal. If it holds ω<ω*<sup>p</sup>* , the dielectric constant is negative, therefore its square root is imaginary pure; this involves the reflection of the incident wave. In the contrary case, the square root of the dielectric constant is real and the

The plasma oscillations are the fluctuations of charge density, which happen to the frequency ω*p* and propagate in the metal. The quanta of such fluctuations inside the volume

There is also a mechanism on the metal surface, characterized by the "surface plasmons". Normally it appears in the wavelength range between the visible and the infrared for the interface air/metal. Localized surface electromagnetic excitations can exist at the surface of nanoparticles and metallic nanostructures: the "localized surface plasmons" (LSP) [11]. The frequency and the intensity of the radiation are very sensitive with respect to the dimension, form and morphology of the nanostructures [12]. Nanoparticles and nanostructures with smaller dimensions with respect to the wavelength of the exciting light are characterized by a wide absorption band, normally in the range of the visible and

The models concerning plasmonics simulate the extinction of the "localized surface plasmons resonance" (LSPR) from nanoparticles and nanoarrays; the application regards the

The Mie theory is a theoretical approach concerning the optical properties of the nanoparticles. When the nanoparticle dimension is smaller than the wavelength of the

calculation of the light absorption and the scattering. The most used models are:

• the discrete-dipole approximation method (DDA) [15]; • the finite-difference time-domain method (FDTD) [16].

incident wave can propagate in the medium with a small attenuation.

The Gans theory extends the Mie theory to spheroidal particles case. In the case of polarization of incident light parallel to the symmetry axis of the spheroid, the cross section of the process is given by:

$$\sigma\_{ext} = \frac{2\pi V \,\mathrm{e}\_{m}^{3/2}}{3\lambda} \sum\_{j} \left[ \frac{(1/P\_{j})^{2} \,\mathrm{e}\_{j}}{(\mathrm{e}\_{r} + \frac{1-P\_{j}}{P\_{j}} \mathrm{e}\_{m})^{2} + \mathrm{e}\_{j}^{2}} \right] \tag{15}$$

with *V* volume of nanoparticle and *Pj* depolarization factors along the three cartesian aces, which hold into account of the anisotropic form of particles.

#### *6.3.3. The discrete-dipole approximation method (DDA)*

With this method the nanoparticles are divided in a cubic array of *N* polarizable dipoles, with polarizabilities *<sup>i</sup>* α determined by the dielectric function of nanoparticles. The induced dipole *<sup>i</sup> P* of every element results *i i loc i*, *P E* =α in an applied plane wave field. The cross section is determined by:

$$\sigma\_{ext} = \frac{4\,\pi k}{\left|\vec{E}\_{inc}\right|^2} \sum\_{j=1}^{N} \text{Im}\left(\vec{E}\_{inc,j} \stackrel{\*}{\cdot} \vec{P}\_j\right) \tag{16}$$

with *k c* =ω (in vacuum). Such method can be applied for the calculation of the absorption, scattering, extinction and other optical properties of nanoparticles of various forms and dimensions. It is considered one of the most important methods for the understanding of the structural characteristics and optical properties of nanomaterials and nanostructures.

#### *6.3.4. The finite-difference time-domain method (FDTD)*

This is a method of numerical calculation in order to resolve the Maxwell equations directly in the time domain. Being a time-domain method, the solutions can cover a wide frequency range with a single process of simulation. It is a versatile and useful technique in applications where the resonance frequencies are not exactly known. A great variety of magnetic and dielectric materials can be modelled in relatively simple way with such method [17].

#### **7. Linear response theory: a new interesting idea**

We start considering a system with an hamiltonian of the form:

$$H = H\_0 + H\_1 \tag{17}$$

Modelling at Nanoscale 11

Analytical calculations permit to write the real part of the complex conductivity σ ω( ) as:

Re ( ) ( ) 1 *<sup>e</sup> KT S e V*

( ) (0) ( ) *i t*

*S dt r r t e* + ∞ α β − ω

The part ..... *<sup>T</sup>* in the integral (24) is the thermal average, and the exponential factor arises

, *d i v r Hr dt*

= = ,

Eq. (23) can be written in a form containing the velocities correlation function instead of the position correlation function. Assuming the high temperature limit ω<<*KT* (valid in such

The integral in Eq. (25) spans the entire *t* axis, so we can perform the complete inverse

*KTV v vt d e e*

< >= ω σ ω

<sup>2</sup> (0) ( ) Re ( ) *i t*

The new introduced key idea is the possibility to perform a complete inversion of Eq. (26) on temporal scale, i.e. considering the entire time axis (-∞, +∞), not the half time axis (0, +∞), as usually considered in literature [9,18]. This idea is viable if we consider the real part of the complex conductivity σ ω( ) . Via contour integration by the residue theorem in Eq. (26), the integral is determined by the poles of Re ( ) σ ω . This leads to an exact formulation and gives a powerful method to describe the velocities correlation function (and consequently the

Another interesting quantity at nanolevel is the mean squared displacement of particles at

2 2 *R t Rt R* ( ) [ ( ) (0)] = −

+ ∞ α β ω

−∞

*<sup>e</sup> dt v v t e*

+ ∞ α β − ω

βα βα ωπ σ ω= ω −

from equilibrium thermal weights for Fermi particles. By considering the identity

2 Re ( ) (0) ( ) <sup>2</sup>

*V KT*

*T*

mean square displacement and the diffusion coefficient) in analytical way.

**8. The other important functions in the nano-bio-context** 

with:

contests), we obtain:

equilibrium, defined as:

Fourier transform of this equation. It gives:

( ) <sup>2</sup>

− ω

*T*

βα −∞ ω = (24)

*i t*

(27)

*T*

βα −∞ σ ω= (25)

βα

<sup>π</sup> (26)

(23)

with *H*1 having small effects with respect to *H*<sup>0</sup> , and of the form:

$$H\_1 = \lambda \, A \, e^{-i \, \alpha t} \, e^{\eta t} \tag{18}$$

being λ a real quantity and η positive, so that in remote past the perturbation is negligible (adiabatic representation: 1 lim 0 *<sup>t</sup>*→−∞ *H* = ). In the case of an electric field of frequency ω we have:

$$H\_1 = e \, \vec{E} \cdot \vec{r} \tag{19}$$

If the electric field is constant in space and its evolution depends on time, it is writable as follows:

$$
\vec{E} = \vec{E}\_0 \, e^{i\alpha t} \tag{20}
$$

The time dependent corresponding current is:

$$
\vec{f}(t) = \sigma(\text{or})\,\vec{E}(t) \tag{21}
$$

The conductivity σ ω( ) is in general a complex function of the frequency ω and can be deduced from linear response theory. Following the standard time-dependent approach, it is possible to find a general formula for the linear response of a dipole moment density *B er V* = in the β direction with the electric field *<sup>E</sup>* in the α direction, where *V* is the volume of the system. This permits to deduce the susceptivity χ( ) ω , which is correlated to σ ω( ) through the relation:

$$1 + 4\,\pi\chi(\alpha) = 1 + 4\,\pi i \frac{\sigma(\alpha)}{\alpha} \tag{22}$$

Analytical calculations permit to write the real part of the complex conductivity σ ω( ) as:

$$\operatorname{Re}\sigma\_{\beta\alpha}(\text{co}) = \frac{e^2 \,\alpha\pi}{V\,\hbar} S\_{\beta\alpha}(\text{co}) \left(1 - e^{-\hbar\alpha\beta KT}\right) \tag{23}$$

with:

10 Plasmonics – Principles and Applications

method [17].

have:

follows:

*B er V* =

σ ω( ) through the relation:

*6.3.4. The finite-difference time-domain method (FDTD)* 

**7. Linear response theory: a new interesting idea** 

We start considering a system with an hamiltonian of the form:

with *H*1 having small effects with respect to *H*<sup>0</sup> , and of the form:

The time dependent corresponding current is:

in the β direction with the electric field *<sup>E</sup>*

1

being λ a real quantity and η positive, so that in remote past the perturbation is negligible (adiabatic representation: 1 lim 0 *<sup>t</sup>*→−∞ *H* = ). In the case of an electric field of frequency ω we

*H eE r* <sup>1</sup> = ⋅

If the electric field is constant in space and its evolution depends on time, it is writable as

0 *i t E Ee* <sup>ω</sup> = 

*<sup>J</sup>* () ( ) ( ) *t Et* =σ ω

The conductivity σ ω( ) is in general a complex function of the frequency ω and can be deduced from linear response theory. Following the standard time-dependent approach, it is possible to find a general formula for the linear response of a dipole moment density

volume of the system. This permits to deduce the susceptivity χ( ) ω , which is correlated to

( ) 14 ()14 *<sup>i</sup>*

+ πχ ω = + π

σ ω

This is a method of numerical calculation in order to resolve the Maxwell equations directly in the time domain. Being a time-domain method, the solutions can cover a wide frequency range with a single process of simulation. It is a versatile and useful technique in applications where the resonance frequencies are not exactly known. A great variety of magnetic and dielectric materials can be modelled in relatively simple way with such

*HH H* 0 1 = + (17)

*it t H Ae e* −ω η =λ (18)

(19)

(20)

(21)

in the α direction, where *V* is the

<sup>ω</sup> (22)

$$S\_{\beta\alpha}(\alpha) = \int\_{-\infty}^{+\infty} dt \left\langle \vec{r}^{\alpha}(0) \,\vec{r}^{\beta}(t) \right\rangle\_{T} e^{-i\alpha t} \tag{24}$$

The part ..... *<sup>T</sup>* in the integral (24) is the thermal average, and the exponential factor arises from equilibrium thermal weights for Fermi particles. By considering the identity

$$v = \frac{d}{dt}r = \frac{i}{\hbar} \left[H\_{\prime\prime}r\right]\_{\prime\prime}$$

Eq. (23) can be written in a form containing the velocities correlation function instead of the position correlation function. Assuming the high temperature limit ω<<*KT* (valid in such contests), we obtain:

$$\operatorname{Re}\sigma\_{\beta\alpha}(\mathbf{u}) = \frac{e^2}{2V\,KT} \int\_{-\infty}^{+\infty} dt \left\langle \bar{v}^{\alpha}(0)\bar{v}^{\beta}(t) \right\rangle\_T e^{-i\alpha t} \tag{25}$$

The integral in Eq. (25) spans the entire *t* axis, so we can perform the complete inverse Fourier transform of this equation. It gives:

$$<\bar{v}^{\alpha}(0)\,\,\bar{v}^{\beta}(t)>\_{T} = \frac{KT}{\pi e^{2}}\int\_{-\infty}^{+\infty} d\alpha \,\,\mathrm{Re}\,\sigma\_{\beta\alpha}(\alpha)\,e^{i\alpha t}\tag{26}$$

The new introduced key idea is the possibility to perform a complete inversion of Eq. (26) on temporal scale, i.e. considering the entire time axis (-∞, +∞), not the half time axis (0, +∞), as usually considered in literature [9,18]. This idea is viable if we consider the real part of the complex conductivity σ ω( ) . Via contour integration by the residue theorem in Eq. (26), the integral is determined by the poles of Re ( ) σ ω . This leads to an exact formulation and gives a powerful method to describe the velocities correlation function (and consequently the mean square displacement and the diffusion coefficient) in analytical way.

#### **8. The other important functions in the nano-bio-context**

Another interesting quantity at nanolevel is the mean squared displacement of particles at equilibrium, defined as:

$$R^2(t) = \left\langle \left[\vec{R}(t) - \vec{R}(0)\right]^2 \right\rangle \tag{27}$$

#### 12 Plasmonics – Principles and Applications

where *R t*( ) is the position vector at time *t*. Through a coordinate transformation relative to the integration region it is possible to rewrite the relation (27) as follows:

$$R^2(t) = 2\int\_0^t dt' \langle t - t' \rangle \left\langle \vec{v}(t') \cdot \vec{v}(0) \right\rangle \tag{28}$$

Modelling at Nanoscale 13

**10. About the Drude-Lorentz-like models: a new promising "plasmon** 

materials, but it holds also for other objects, like ions [19].

quantum behaviour are the weights *<sup>i</sup>*

*m*

decaying times of each mode respectively).

Recently it has been published a new formulation of the Drude-Lorentz model [1,20], based on linear response theory and resonant plasmonic mode; it is able to accommodate some of the observed departures and gives a detailed description of the dynamic response of the carriers at nano-level. One of the peculiarities of this new model is the inversion of relation (26), extending the integration on time to the entire time axis (-∞, +∞); this procedure is not trivial and it was introduced for the first time in the nano-bio-context. One of the main advantages is the possibility to obtain analytical relations for the description of the dynamic behaviour of nanosystems. Such relations are functions of parameters, experimentally obtainable through Time-resolved techniques. An important consequence of the application of such formulation turns out to be the mathematical justification of the unexpected experimental results of high initial mobility of charge carriers in devices based on semiconductor nanomaterials, as the dye sensitized solar cells (DSSCs). The increase of the diffusion coefficient implies a raise of the efficiency of such electro-chemical cells. The new formulation is able to describe very adequately the properties of transport of nano-bio-

Many experimental data have indicated that plasmon models describe nanostructured systems in particularly effective way. At quantum level, the key factors incorporating the

> <sup>2</sup> 4 *<sup>i</sup> p i N e <sup>f</sup> m* π

<sup>1</sup> (0) ( ) exp cos sin

*KT tt t v vt <sup>f</sup>*

( ) ( ) (1 ) (1 ) <sup>1</sup> (0) ( ) 1 exp 1 exp 2 22

*KT <sup>f</sup> t t v vt*

 α α ⋅ = <sup>−</sup> <sup>−</sup> τ τα τ (36)

with: 2 2 4 1 *iR i i* α = τω− ; (37)

 +α −α ⋅ = +α − − −α − α ττ (38)

with: 2 2 1 4 *i I i i* α= −τω . (39)

(*K* is the Boltzmann's constant, *T* the temperature of the system, ω*i* and *<sup>i</sup>* τ frequencies and

The velocities correlation function in quantum case assumes the form:

*i*

*m*

where *N* is the carrier density, *m* and *e* respectively the mass and the charge of the electron.

*f* , related to plasma frequencies by:

22 2

*i i i iR i*

*i I i I i i I i I i i I i i*

ω = (35)

*i R i R*

2

**model"** 

The mean squared displacement can therefore be evaluated through the velocities correlation function. Through Eq. (26) it is possible to deduce also the diffusion coefficient *D*:

$$D(t) = \frac{d\mathbb{R}^2(t)}{2\,dt} = \int\_0^t dt' \left\langle \vec{v}(t') \cdot \vec{v}(0) \right\rangle \tag{29}$$

The three relations (26), (28), (29) are fundamental in deducing the most important characteristics concerning the transport phenomena.

#### **9. About the complex conductivity** σ ω( )

Let us consider a particle in a region in which there is an electric time-oscillating field, direct along *z*-axis, with an elastic-type and a friction-type force acting on system; the dynamic equation of the particle can be written as:

$$m\ddot{z} = -k\,z - \lambda\,\dot{z} + e\,E\_o\,e^{-i\,\alpha t} \tag{30}$$

where *m* is the mass of the particle, *k* the strenght constant of the oscillators, λ= τ *m* and *τ* is the relaxation time. We consider solutions of the form:

$$z(t) = z\_0 e^{-i\alpha t} \tag{31}$$

The current density in the field direction is:

$$j\_z = \imath e \,\dot{z} = \imath e \, z\_0 \left( -i \,\text{o} \right) e^{i \,\text{o} \, t} = \mathbf{c} \, E \tag{32}$$

From Eq. (32), via analytical calculations, it is possible to obtain the frequency dependent complex conductivity:

$$\sigma = \frac{i\,\alpha m e^2}{m(\,\alpha^2 - \alpha\_0^2 + i\,\lambda\,\alpha)}\tag{33}$$

The real part of the complex conducivity (33) results [1,19]:

$$\text{Re}\,\mathbf{\sigma}(\text{co}) = \frac{\pi m e^2}{m} \frac{\text{co}^2/\text{r}^2}{\left(\text{co}\_0^2 - \text{o}^2\right)^2 + \text{o}^2/\text{r}^2} \tag{34}$$

### **10. About the Drude-Lorentz-like models: a new promising "plasmon model"**

Recently it has been published a new formulation of the Drude-Lorentz model [1,20], based on linear response theory and resonant plasmonic mode; it is able to accommodate some of the observed departures and gives a detailed description of the dynamic response of the carriers at nano-level. One of the peculiarities of this new model is the inversion of relation (26), extending the integration on time to the entire time axis (-∞, +∞); this procedure is not trivial and it was introduced for the first time in the nano-bio-context. One of the main advantages is the possibility to obtain analytical relations for the description of the dynamic behaviour of nanosystems. Such relations are functions of parameters, experimentally obtainable through Time-resolved techniques. An important consequence of the application of such formulation turns out to be the mathematical justification of the unexpected experimental results of high initial mobility of charge carriers in devices based on semiconductor nanomaterials, as the dye sensitized solar cells (DSSCs). The increase of the diffusion coefficient implies a raise of the efficiency of such electro-chemical cells. The new formulation is able to describe very adequately the properties of transport of nano-biomaterials, but it holds also for other objects, like ions [19].

Many experimental data have indicated that plasmon models describe nanostructured systems in particularly effective way. At quantum level, the key factors incorporating the quantum behaviour are the weights *<sup>i</sup> f* , related to plasma frequencies by:

$$\left|\alpha\_{p\_i}\right|^2 = \frac{4\pi N e^2}{m} f\_i \tag{35}$$

where *N* is the carrier density, *m* and *e* respectively the mass and the charge of the electron.

The velocities correlation function in quantum case assumes the form:

$$
\left\langle \vec{v}(0) \cdot \vec{v}(t) \right\rangle = \left(\frac{KT}{m}\right) \sum\_{i} \left( f\_i \exp\left(-\frac{t}{2\tau\_i}\right) \left[ \cos\left(\frac{\alpha\_{iR}}{2} \frac{t}{\tau\_i}\right) - \frac{1}{\alpha\_{iR}} \sin\left(\frac{\alpha\_{iR}}{2} \frac{t}{\tau\_i}\right) \right] \right) \tag{36}
$$

12 Plasmonics – Principles and Applications

*D*:

where *R t*( ) is the position vector at time *t*. Through a coordinate transformation relative to

( ) 2 ' ' ( ') (0)

The mean squared displacement can therefore be evaluated through the velocities correlation function. Through Eq. (26) it is possible to deduce also the diffusion coefficient

> 0 ( ) ( ) ' ( ') (0) <sup>2</sup> *<sup>t</sup> dR t D t dt v t v*

The three relations (26), (28), (29) are fundamental in deducing the most important

Let us consider a particle in a region in which there is an electric time-oscillating field, direct along *z*-axis, with an elastic-type and a friction-type force acting on system; the dynamic

*i t mz kz z eE e <sup>o</sup>*

where *m* is the mass of the particle, *k* the strenght constant of the oscillators, λ= τ *m* and *τ* is

<sup>0</sup> ( ) *i t*

2

( )

0

τ ωτ

<sup>2</sup> 2 2 22

ω −ω +ω τ

2 22

From Eq. (32), via analytical calculations, it is possible to obtain the frequency dependent

2 2 <sup>0</sup> ( ) *i ne m i* ω

*R t dt t t v t v* =− ⋅ (28)

*dt* == ⋅ (29)

− ω =− −λ + (30)

<sup>0</sup> ( ) *i t zt z e*− ω <sup>=</sup> (31)

ω −ω + λω (33)

(34)

*j nez nez i e E* <sup>ω</sup> = = − ω =σ (32)

the integration region it is possible to rewrite the relation (27) as follows:

( ) <sup>2</sup> 0

2

characteristics concerning the transport phenomena.

the relaxation time. We consider solutions of the form:

*z*

The real part of the complex conducivity (33) results [1,19]:

σ=

Re ( ) *ne*

σω=

*m*

**9. About the complex conductivity** σ ω( )

equation of the particle can be written as:

The current density in the field direction is:

complex conductivity:

*t*

$$\alpha\_{iR} = \sqrt{4\,\tau\_i^2 \,\mathrm{o}\_i^2 - 1} \;:\tag{37}$$

$$\left\langle \vec{v}(0) \cdot \vec{v}(t) \right\rangle = \frac{1}{2} \left( \frac{KT}{m} \right) \sum\_{i} \left[ \left( \frac{f\_i}{\alpha\_{iI}} \right) \left[ \left( 1 + \alpha\_{iI} \right) \exp \left( -\frac{\left(1 + \alpha\_{iI}\right)}{2} \frac{t}{\tau\_i} \right) - \left( 1 - \alpha\_{iI} \right) \exp \left( -\frac{\left(1 - \alpha\_{iI}\right)}{2} \frac{t}{\tau\_i} \right) \right] \right] (38)$$

$$\alpha\_{iI} = \sqrt{1 - 4\,\tau\_i^2 \,\mathrm{o}\_i^2} \,\, . \tag{39}$$

(*K* is the Boltzmann's constant, *T* the temperature of the system, ω*i* and *<sup>i</sup>* τ frequencies and decaying times of each mode respectively).

Integrating these expressions through Eqs. (28) and (29), we obtain the mean squared displacement <sup>2</sup> *R* and the diffusion coefficient *D* [20,21].

Modelling at Nanoscale 15

This means that diffusion occurs after sufficient time has elapsed, so that scattering events

<sup>0</sup> ω= ⋅ 1.12 10 *Hz* dot-dashed line; <sup>11</sup>

<sup>0</sup> ω= ⋅ 2.24 10 *Hz*

become significant, while at smaller times the motion is essentially ballistic.

**Figure 2.** *R*2 vs *x t* = τ for 2 values of τ ( <sup>11</sup>

dashed line) for TiO2 (*m*=6*m*e, *T*=300K). Saturation values occur at sufficiently large *t*.

mobility and to the existence of a strictly insulating state for t→∞

Figure 3, with the scattering time in the range 10-13-10-14*s*.

The Drude behaviour is contrasted by the plasmon behaviour. It has reported representative cases at large times *x t* = τ>>1 (Figure 2) for TiO2 nanoparticle films with small frequencies ω<sup>0</sup> , i.e. close to the Drude case. It is observable that all of the curves reach a plateau value at sufficiently long times and that the slope within a given time interval increases with τ . As consequences of this behaviour, the plateau of <sup>2</sup> *R* may become larger than the size of the nanoparticles composing the films, depending on the parameters ω0 and τ , indicating that carriers created at time 0 *t* = will have enhanced mobility in the nanoporous films at small times, on the order of few τ , in contrast with a commonly expected low mobility in the disordered TiO2 network. Secondly, the time derivative of <sup>2</sup> *R* decreases to zero at very large times, indicating localization by scattering and therefore leading to decreasing

Few cases corresponding to short time behaviour and large frequencies are showed in

### **11. Some interesting results of the new model**

The model has demonstrated at today high generality and good accordance with experiments. It has been applied in relation to the most commonly used and studied materials at nano-level, i.e. zinc oxide (ZnO), titanium dioxide (TiO2), gallium arsenide (GaAs), silicon (Si) and carbon nanotubes (CN), also at nano-bio-sensoristic level [19,22].

In Figure 1 it is showed <sup>2</sup> *R* for doped Silicon. For this systems the conductivity is the contribution of a Drude-Lorentz term and a Drude term. At large times the Drude-Lorentz term leads to an <sup>2</sup> *R* approaching a constant value (Figure 2), while the Drude term alone is the dominant one. Therefore, for sufficiently large times, only the Drude term survives.

**Figure 1.** *R*2 vs *x t* = τ for some representative values of τ , typical of doped Silicon [23] ( <sup>0</sup> ω =0 , *T*=300K) (Drude model). A complete description of *R*2 for Silicon requires the evaluation of the contribution of the Drude-Lorentz part.

We can observe that the linear relation at large times becomes quadratic at smaller times. The cross-over between the two regimes occurs at times comparable to the scattering time. This means that diffusion occurs after sufficient time has elapsed, so that scattering events become significant, while at smaller times the motion is essentially ballistic.

14 Plasmonics – Principles and Applications

displacement <sup>2</sup> *R* and the diffusion coefficient *D* [20,21].

**11. Some interesting results of the new model** 

Integrating these expressions through Eqs. (28) and (29), we obtain the mean squared

The model has demonstrated at today high generality and good accordance with experiments. It has been applied in relation to the most commonly used and studied materials at nano-level, i.e. zinc oxide (ZnO), titanium dioxide (TiO2), gallium arsenide (GaAs), silicon (Si) and carbon nanotubes (CN), also at nano-bio-sensoristic level [19,22].

In Figure 1 it is showed <sup>2</sup> *R* for doped Silicon. For this systems the conductivity is the contribution of a Drude-Lorentz term and a Drude term. At large times the Drude-Lorentz term leads to an <sup>2</sup> *R* approaching a constant value (Figure 2), while the Drude term alone is the dominant one. Therefore, for sufficiently large times, only the Drude term survives.

**Figure 1.** *R*2 vs *x t* = τ for some representative values of τ , typical of doped Silicon [23] ( <sup>0</sup> ω =0 , *T*=300K) (Drude model). A complete description of *R*2 for Silicon requires the evaluation of the

We can observe that the linear relation at large times becomes quadratic at smaller times. The cross-over between the two regimes occurs at times comparable to the scattering time.

contribution of the Drude-Lorentz part.

**Figure 2.** *R*2 vs *x t* = τ for 2 values of τ ( <sup>11</sup> <sup>0</sup> ω= ⋅ 1.12 10 *Hz* dot-dashed line; <sup>11</sup> <sup>0</sup> ω= ⋅ 2.24 10 *Hz* dashed line) for TiO2 (*m*=6*m*e, *T*=300K). Saturation values occur at sufficiently large *t*.

The Drude behaviour is contrasted by the plasmon behaviour. It has reported representative cases at large times *x t* = τ>>1 (Figure 2) for TiO2 nanoparticle films with small frequencies ω<sup>0</sup> , i.e. close to the Drude case. It is observable that all of the curves reach a plateau value at sufficiently long times and that the slope within a given time interval increases with τ . As consequences of this behaviour, the plateau of <sup>2</sup> *R* may become larger than the size of the nanoparticles composing the films, depending on the parameters ω0 and τ , indicating that carriers created at time 0 *t* = will have enhanced mobility in the nanoporous films at small times, on the order of few τ , in contrast with a commonly expected low mobility in the disordered TiO2 network. Secondly, the time derivative of <sup>2</sup> *R* decreases to zero at very large times, indicating localization by scattering and therefore leading to decreasing mobility and to the existence of a strictly insulating state for t→∞

Few cases corresponding to short time behaviour and large frequencies are showed in Figure 3, with the scattering time in the range 10-13-10-14*s*.

Modelling at Nanoscale 17

**Figure 4.** Velocities correlation function vs *x t* = τ for some values of *<sup>I</sup>* α (*m*=6*m*e, *T*=300K).

[27].

typical for bulk GaAs at room temperature.

In systems with charge localization, where the carrier mean free path is comparable to the characteristic dimension of the nanoparticles, the conductivity response becomes more complicated. The response at low frequencies is characterized by an increasing real part of the conductivity and by a negative imaginary part. Deviations from the Drude model become strong in nanostructured materials, such as photoexcited TiO2 nanoparticles, ZnO films, InP nanoparticles [25], semiconducting polymer molecules [26], and carbon nanotubes

In isolated GaAs nanowires the electronic response exhibits a pronounced surface plasmon mode, that forms within 300 fs, before decaying within 10 ps as a result of charge trapping at the nanowire surface. The conductivity in this case was fitted by using the Drude model for a plasmon and the mobility found to be remarkably high, being roughly one-third of that

The Smith model with 1*c* = −1 is obtained as a limit of the new introduced plasmon model when 0 *<sup>I</sup>* α → . On performing this limit in Eq. (38) one finds the expression obtained by Smith with a scattering time twice the plasmon scattering time τ . The situation 0 *<sup>I</sup>* α → corresponds to 1 / 2 *<sup>o</sup>* τ≈ ω . From the other hand, both Smith's and the new model reduce to the Drude model in the limit 0 ω →*<sup>o</sup>* . So, although the two models are analytically different, their predictions are expected to be quite similar. From Figure 4 it is possible to note that, for

**Figure 3.** *R*2 vs *x t* = τ for 3 values of τ ( <sup>13</sup> <sup>0</sup> ω= ⋅ 0.5 10 *Hz* solid line; <sup>13</sup> <sup>0</sup> ω =10 *Hz* dot-dashed line; 14 <sup>0</sup> ω= ⋅ 0.5 10 *Hz* dashed line) for TiO2 (*m*=6*m*e, *T*=300K).

In relation to the behaviour of the velocities correlation function, excluding the Drude model ( <sup>0</sup> ω = 0 ), characterized by simple exponential decreasing functions with time, the correlation function of velocities corresponds to either a damped oscillatory behaviour (Eq. (36)), or to a superposition of two exponentials (Eq. (38)), with quite different decay times depending on the value of *<sup>I</sup>* α (Figures 4, 5).

At 1 α ≅ , i.e. close to Drude behaviour, proper case of some systems [3], the current results the superposition of a short ≈ τ and long / (1 )*<sup>I</sup>* τ −α time decay modes. The current is a damped oscillating function of time (Eq. (36)) when 0 ω τ>>1 , and will have double exponential behaviour for 0 ω τ<<1 (Eq. (38)).

From this analysis, it is viable the possibility that the previous results can give an explanation of the ultra-short times and high mobilities, with which the charges spread in mesoporous systems, of large interest in photocatalitic and photovoltaic systems [24]; the relative short times (few τ ), with which charges can reach much larger distances than typical dimensions of nanoparticles, indicate easy charges diffusion inside the nanoparticles. The unexplained fact, experimentally found, of ultrashort injection of charge carriers (particularly in Grätzel's cells) can be related to this phenomenon.

16 Plasmonics – Principles and Applications

**Figure 3.** *R*2 vs *x t* = τ for 3 values of τ ( <sup>13</sup>

depending on the value of *<sup>I</sup>* α (Figures 4, 5).

exponential behaviour for 0 ω τ<<1 (Eq. (38)).

(particularly in Grätzel's cells) can be related to this phenomenon.

<sup>0</sup> ω= ⋅ 0.5 10 *Hz* dashed line) for TiO2 (*m*=6*m*e, *T*=300K).

14

<sup>0</sup> ω= ⋅ 0.5 10 *Hz* solid line; <sup>13</sup>

In relation to the behaviour of the velocities correlation function, excluding the Drude model ( <sup>0</sup> ω = 0 ), characterized by simple exponential decreasing functions with time, the correlation function of velocities corresponds to either a damped oscillatory behaviour (Eq. (36)), or to a superposition of two exponentials (Eq. (38)), with quite different decay times

At 1 α ≅ , i.e. close to Drude behaviour, proper case of some systems [3], the current results the superposition of a short ≈ τ and long / (1 )*<sup>I</sup>* τ −α time decay modes. The current is a damped oscillating function of time (Eq. (36)) when 0 ω τ>>1 , and will have double

From this analysis, it is viable the possibility that the previous results can give an explanation of the ultra-short times and high mobilities, with which the charges spread in mesoporous systems, of large interest in photocatalitic and photovoltaic systems [24]; the relative short times (few τ ), with which charges can reach much larger distances than typical dimensions of nanoparticles, indicate easy charges diffusion inside the nanoparticles. The unexplained fact, experimentally found, of ultrashort injection of charge carriers

<sup>0</sup> ω =10 *Hz* dot-dashed line;

**Figure 4.** Velocities correlation function vs *x t* = τ for some values of *<sup>I</sup>* α (*m*=6*m*e, *T*=300K).

In systems with charge localization, where the carrier mean free path is comparable to the characteristic dimension of the nanoparticles, the conductivity response becomes more complicated. The response at low frequencies is characterized by an increasing real part of the conductivity and by a negative imaginary part. Deviations from the Drude model become strong in nanostructured materials, such as photoexcited TiO2 nanoparticles, ZnO films, InP nanoparticles [25], semiconducting polymer molecules [26], and carbon nanotubes [27].

In isolated GaAs nanowires the electronic response exhibits a pronounced surface plasmon mode, that forms within 300 fs, before decaying within 10 ps as a result of charge trapping at the nanowire surface. The conductivity in this case was fitted by using the Drude model for a plasmon and the mobility found to be remarkably high, being roughly one-third of that typical for bulk GaAs at room temperature.

The Smith model with 1*c* = −1 is obtained as a limit of the new introduced plasmon model when 0 *<sup>I</sup>* α → . On performing this limit in Eq. (38) one finds the expression obtained by Smith with a scattering time twice the plasmon scattering time τ . The situation 0 *<sup>I</sup>* α → corresponds to 1 / 2 *<sup>o</sup>* τ≈ ω . From the other hand, both Smith's and the new model reduce to the Drude model in the limit 0 ω →*<sup>o</sup>* . So, although the two models are analytically different, their predictions are expected to be quite similar. From Figure 4 it is possible to note that, for

#### 18 Plasmonics – Principles and Applications

0 *<sup>I</sup>* α → , the lower curve is of the form of the Smith curve. The backscattering mechanism invoked by Smith arises in a natural way in this new model, without further assumptions on successive scattering events.

Modelling at Nanoscale 19

<sup>0</sup> ω≅ ⋅ 1.5 10 *Hz* ,

<sup>0</sup> ω≅ ⋅ 2.5 10 *Hz* (for

from which we deduce 0 ω τ≅0.2 . For τ of the order of 10-13*s*, we obtain <sup>12</sup>

procedure for injection times of films (500 nm grains) leads to <sup>12</sup>

<sup>13</sup> 10 *s* <sup>−</sup> τ ≅ ).

<sup>0</sup> ω τ= ÷ 0.3 0.49 .

model.

interval [0.1-1] and

τ

of the new plasmon model.

which is of the correct order of magnitude for the resonance in the infrared. The same

In GaAs nanowires [27,28] a resonance at 0.3 0.5 *<sup>o</sup>* ω= − *THz* has been suggested to explain the frequency-dependent conductivity obtained from THz experiments. A long characteristic obtained time <sup>12</sup> 1.1 5 10 *<sup>c</sup> <sup>s</sup>* <sup>−</sup> τ≈ ÷⋅ is in accordance with the new model, corresponding to <sup>12</sup> / (1 ) 5.10 *c I <sup>s</sup>* <sup>−</sup> τ =τ −α = and leading to 0.1 0.8 *<sup>I</sup>* α= ÷ for <sup>12</sup> <sup>10</sup> *<sup>s</sup>* <sup>−</sup> τ = and

Similar considerations can be applied to single walled carbon nanotubes (SWCN), where the resonance state is at 0.5 *<sup>o</sup>* ω = *THz* with Drude-Lorentz scattering time τ ≈ 10 -13 s and to sensitized TiO2 nanostructured films. In the latter the conductivity of electrons injected from the excited state of the dye molecule into the conduction band of TiO2 is initially more Drude-like, and then evolves into a conduction dominated by strong backscattering as the carriers equilibrate with the lattice. These results conform to the backscattering mechanism

In Figure 6 it is presented the fitting of data of GaAs photoconductivity [27,28] with the new

**Figure 6.** Behaviour of the diffusion coefficient *D* for *<sup>I</sup>* α (Eq. (39) at classical level [1]) varying in the

= 10-12 *s*. Note that *D*→0 for *t*→∞, indicating absence of diffusion at very long times.

**Figure 5.** Velocities correlation function vs *x t* = τ for two values of *<sup>R</sup>* α (*m*=6*m*e, *T*=300K). Evident exponentially damped oscillations are displayed in this case.

Time evolution of charges can be systematically studied by means of time resolved techniques. The THz technique allows a detailed investigation of very short time behaviour, such as the photoinjection, as well as longer time behaviour such as thermalized motion. The studies have indicated a common mechanism of the short time domain in nanostructures, where carriers are close to Drude behaviour with a rather large diffusion coefficient, followed by a range with time decreasing mobility. This latter stage is characterized by decay of the response as the superposition of a short time and a longer time exponentials.

The time response of ZnO films, nanowires, and nanoparticles to near-UV photoexcitation has been investigated in THz experiments. Films and nanoparticles show ultrafast injection, but with the addition of a second slower component. For ZnO nanoparticles, double exponential decay indicates characteristic times 1τ ≈ 94 *ps* and 2 τ ≈ 2.4*ns* . From Eq. (38) it is possible to deduce a ratio of two different relaxation times 1τ and 2 τ ; it is 1 2 / (1 ) / (1 ) *I I* τ τ = −α +α . Therefore, on using the experimental values, we find 0.92 *<sup>I</sup>* α = , from which we deduce 0 ω τ≅0.2 . For τ of the order of 10-13*s*, we obtain <sup>12</sup> <sup>0</sup> ω≅ ⋅ 1.5 10 *Hz* , which is of the correct order of magnitude for the resonance in the infrared. The same procedure for injection times of films (500 nm grains) leads to <sup>12</sup> <sup>0</sup> ω≅ ⋅ 2.5 10 *Hz* (for <sup>13</sup> 10 *s* <sup>−</sup> τ ≅ ).

18 Plasmonics – Principles and Applications

successive scattering events.

0 *<sup>I</sup>* α → , the lower curve is of the form of the Smith curve. The backscattering mechanism invoked by Smith arises in a natural way in this new model, without further assumptions on

**Figure 5.** Velocities correlation function vs *x t* = τ for two values of *<sup>R</sup>* α (*m*=6*m*e, *T*=300K). Evident

Time evolution of charges can be systematically studied by means of time resolved techniques. The THz technique allows a detailed investigation of very short time behaviour, such as the photoinjection, as well as longer time behaviour such as thermalized motion. The studies have indicated a common mechanism of the short time domain in nanostructures, where carriers are close to Drude behaviour with a rather large diffusion coefficient, followed by a range with time decreasing mobility. This latter stage is characterized by decay of the response as the superposition of a short time and a longer time

The time response of ZnO films, nanowires, and nanoparticles to near-UV photoexcitation has been investigated in THz experiments. Films and nanoparticles show ultrafast injection, but with the addition of a second slower component. For ZnO nanoparticles, double exponential decay indicates characteristic times 1τ ≈ 94 *ps* and 2 τ ≈ 2.4*ns* . From Eq. (38) it is possible to deduce a ratio of two different relaxation times 1τ and 2 τ ; it is 1 2 / (1 ) / (1 ) *I I* τ τ = −α +α . Therefore, on using the experimental values, we find 0.92 *<sup>I</sup>* α = ,

exponentially damped oscillations are displayed in this case.

exponentials.

In GaAs nanowires [27,28] a resonance at 0.3 0.5 *<sup>o</sup>* ω= − *THz* has been suggested to explain the frequency-dependent conductivity obtained from THz experiments. A long characteristic obtained time <sup>12</sup> 1.1 5 10 *<sup>c</sup> <sup>s</sup>* <sup>−</sup> τ≈ ÷⋅ is in accordance with the new model, corresponding to <sup>12</sup> / (1 ) 5.10 *c I <sup>s</sup>* <sup>−</sup> τ =τ −α = and leading to 0.1 0.8 *<sup>I</sup>* α= ÷ for <sup>12</sup> <sup>10</sup> *<sup>s</sup>* <sup>−</sup> τ = and <sup>0</sup> ω τ= ÷ 0.3 0.49 .

Similar considerations can be applied to single walled carbon nanotubes (SWCN), where the resonance state is at 0.5 *<sup>o</sup>* ω = *THz* with Drude-Lorentz scattering time τ ≈ 10 -13 s and to sensitized TiO2 nanostructured films. In the latter the conductivity of electrons injected from the excited state of the dye molecule into the conduction band of TiO2 is initially more Drude-like, and then evolves into a conduction dominated by strong backscattering as the carriers equilibrate with the lattice. These results conform to the backscattering mechanism of the new plasmon model.

In Figure 6 it is presented the fitting of data of GaAs photoconductivity [27,28] with the new model.

**Figure 6.** Behaviour of the diffusion coefficient *D* for *<sup>I</sup>* α (Eq. (39) at classical level [1]) varying in the interval [0.1-1] and τ= 10-12 *s*. Note that *D*→0 for *t*→∞, indicating absence of diffusion at very long times.

In this figure the upper curve is the Drude result corresponding to 0 ω = 0 . It is the highest value of the diffusion coefficient, as ω0 varies from zero to its maximum value. The dots in Figure 6 are experimental data derived from THz spectroscopy, the error bars on dots represent the distribution of the experimental data and the lines the predictions of the new plasmon model. *D* tends to a vanishing value as *t*→∞, the larger ω0 the faster this vanishing occurs. The results can be interpreted as follows: at early time, of the order of τ, the system behaves as Drude-like, irrespective of ω<sup>0</sup> , with carriers assuming large mobility values (the numerical evaluation gives *D* ≈ 7.5 cm2/s for the case of Figure 6); at increased times, charges become progressively localized as a result of scattering, significantly in agreement with the conclusions drawn by THz spectroscopy [29,30].

Modelling at Nanoscale 21

[8] Han JG, Wan F, Zhu ZY, Liao Y, Ji T, Ge M, Zhang ZY. Shift in low-frequency vibrational spectra of transition-metal zirconium compounds. Applied Physics Letters 2005; 87(17):

[9] Smith NV. Classical generalization of the Drude formula for the optical conductivity.

[10] Shipway AN, Katz E, Willner I. Nanoparticle arrays on surfaces for electronic, optical

[11] Hutter E, Fendler JH. Exploitation of Localized Surface Resonance. Advanced Materials

[12] Willets KA, Van Duyne RP. Localized Surface Plasmon Resonance Spectroscopy and

[13] Stuart DA, Haes AJ, Yonzon CR, Hicks EM, Van Duyne RP. Biological applications of localised surface plasmonic Phenomenae. IEE Proceedings-Nanobiotechnology 2005;

[14] Gulati A, Liao H, Hafner JH. Monitoring Gold Nanorod Synthesis by Localized Surface Plasmon Resonance. The Journal of Physical Chemistry B 2006; 110(45): 22323–22327. [15] Zhang ZY, Zhao Y-P. Optical properties of helical Ag nanostructures calculated by discrete dipole approximation method. Applied Physics Letters 2007; 90: 221501-221503. [16] Sherry LJ, Chang S-H, Schatz GC, Van Duyne RP. Localized Surface Plasmon Resonance Spectroscopy of Single Silver Nanocubes. Nano Letters. 2005; 5(10): 2034–

[17] Zhang W, Yue Z, Wang C, Yang S, Niu W, Liu G. Theoretical models, Fabrications and Applications of Localized Surface Plasmon Resonance sensors. 978-1-4244-4964- 4/10/\$25.00©2010 IEEE 2010; pp. 4; Project supported by the National Natural Science

[18] Sassella A, Borghesi A, Pivac B, Pavesi L. Characterization of porous silicon by microscopic Fourier transform infrared spectroscopy; 266-267. 9th International

[19] Di Sia P. Classical and quantum transport processes in nano-bio-structures: a new theoretical model and applications. PhD Thesis. Verona University – Italy; 2011. [20] Di Sia P. An Analytical Transport Model for Nanomaterials: The Quantum Version.

[21] Di Sia P. Oscillating velocity and enhanced diffusivity of nanosystems from a new

[22] Di Sia P. New theoretical results for high diffusive nanosensors based on ZnO oxides.

[23] Pirozhenko I, Lambrecht A. Influence of slab thickness on the Casimir force. Physical

[24] Grätzel M. Solar Energy Conversion by Dye-Sensitized Photovoltaic Cells. Inorganic

[25] Nienhuys H-K, Sundström V. Influence of plasmons on terahertz conductivity

Journal of Computational and Theoretical Nanoscience 2012; 9: 31-34.

quantum transport model. Journal of Nano Research 2011; 16: 49-54.

measurements. Applied Physics Letters 2005; 87: 012101-012103.

172107-172109.

152(1): 13-32.

2038.

2004; 16(19): 1685-1706.

Physical Review B 2001; 64(15): 155106-155111.

and sensoric applications. ChemPhysChem 2000; 1: 18-52.

Foundation of China (Grant N. 60574091, 60871028).

Conference on Fourier transform spectroscopy; 1994.

Sensors and Transducers Journal 2010; 122(1): 1-8.

Review A 2008; 77: 013811-013818.

Chemistry 2005; 44(20): 6841-6851.

Sensing. Annual Review of Physical Chemistry 2007; 58: 267-297.

### **12. Conclusions**

It is hoped that this new plasmon model will be useful to describe experimental data as an alternative to other generalizations of the Drude model. The principal findings of it are the reversal of the current in small systems like nanoparticles and the time dependence of the transport parameter *D* , which describes the dynamics of the system at short and long time, with increased mobilities at very short times followed by localization at longer times. This unusual behaviour, with respect to the predicted Drude-like, converges as a whole in the peculiar frequency dependence of the optical conductivity. The strength of the new model consist of its ability to accommodate this behaviour and include previous models, like the Smith model [19].

### **Author details**

Paolo Di Sia *Free University of Bozen-Bolzano, Bruneck-Brunico (BK), Italy* 

### **13. References**


[8] Han JG, Wan F, Zhu ZY, Liao Y, Ji T, Ge M, Zhang ZY. Shift in low-frequency vibrational spectra of transition-metal zirconium compounds. Applied Physics Letters 2005; 87(17): 172107-172109.

20 Plasmonics – Principles and Applications

**12. Conclusions** 

Smith model [19].

**Author details** 

**13. References** 

Paolo Di Sia

conclusions drawn by THz spectroscopy [29,30].

*Free University of Bozen-Bolzano, Bruneck-Brunico (BK), Italy* 

and Theoretical Nanoscience 2011; 8: 84-89.

Nanoscience. John Wiley & Sons; 2006.

Science and Technology 2008; 1(1): 1-8.

In this figure the upper curve is the Drude result corresponding to 0 ω = 0 . It is the highest value of the diffusion coefficient, as ω0 varies from zero to its maximum value. The dots in Figure 6 are experimental data derived from THz spectroscopy, the error bars on dots represent the distribution of the experimental data and the lines the predictions of the new plasmon model. *D* tends to a vanishing value as *t*→∞, the larger ω0 the faster this vanishing

behaves as Drude-like, irrespective of ω<sup>0</sup> , with carriers assuming large mobility values (the numerical evaluation gives *D* ≈ 7.5 cm2/s for the case of Figure 6); at increased times, charges become progressively localized as a result of scattering, significantly in agreement with the

It is hoped that this new plasmon model will be useful to describe experimental data as an alternative to other generalizations of the Drude model. The principal findings of it are the reversal of the current in small systems like nanoparticles and the time dependence of the transport parameter *D* , which describes the dynamics of the system at short and long time, with increased mobilities at very short times followed by localization at longer times. This unusual behaviour, with respect to the predicted Drude-like, converges as a whole in the peculiar frequency dependence of the optical conductivity. The strength of the new model consist of its ability to accommodate this behaviour and include previous models, like the

[1] Di Sia P. An Analytical Transport Model for Nanomaterials. Journal of Computational

[2] Wolf EL. Nanophysics and Nanotechnology: An Introduction to Modern Concepts in

[3] Schmuttenmaer CA. Using Terahertz Spectroscopy to Study Nanomaterials. Terahertz

[6] Levy O, Stroud D. Maxwell Garnett theory for mixtures of anisotropic inclusions: Application to conducting polymers. Physical Review B 1997; 56(13): 8035-8046. [7] Han J, Zhang W, Chen W, Ray S, Zhang J, He M, Azad AK, Zhu Z. Terahertz Dielectric Properties and Low-Frequency Phonon Resonances of ZnO Nanostructures. The

[4] Ziman M. Principles of the Theory of Solids. Cambridge University Press; 1979.

[5] Kittel C. Introduction to Solid State Physics. Wiley New York; 1995.

Journal of Physical Chemistry C 2007; 111(35): 13000-13006.

τ

, the system

occurs. The results can be interpreted as follows: at early time, of the order of

	- [26] Hendry E, Koeberg M, Schins JM, Nienhuys HK, Sundström V, Siebbeles LDA, Bonn M. Interchain effects in the ultrafast photophysics of a semiconducting polymer: THz timedomain spectroscopy of thin films and isolated chains in solution. Physical Review B 2005; 71: 125201-125210.

**Chapter 2** 

© 2012 Vandenbosch, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

and reproduction in any medium, provided the original work is properly cited.

© 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

**Computational Electromagnetics in Plasmonics** 

In electromagnetics, numerical techniques have been essential in the development of new technology in the last two decades. The rapidly growing computer capacity and calculation speeds make accurate solutions of very complex problems feasible. This has been especially true in the design of microwave and millimeter wave components and antennas. Whereas 30 years ago, the design of an antenna was based on simple analytical models, or trial and error strategies, nowadays, simulations seem to be as crucial to the

The situation is quite different in plasmonics. Plasmonics is a quite novel research field and the application of computational electromagnetics in plasmonics can be categorized as "very recent". There are many challenges that still need to be faced and "missing links" that have

The plasmonic structures targeted are structures in the order of magnitude of a wavelength at plasmonic frequencies, i.e. at near IR and optical frequencies, and beyond. Although this frequency range is totally different from the traditional range where computational tools have been developed, i.e. the microwave range, in most cases no special modeling techniques have to be used. By far most plasmonic topologies reported in literature have been analyzed / designed with the well-known numerical techniques implemented within in-house developed or commercial software packages. This means that in this chapter the major numerical techniques can be overviewed in a general sense, referring to standard literature. These techniques will not be derived or explained here in full detail. Instead, this chapter focuses on those aspects that come into the picture when the structure is plasmonic. After the section on techniques and tools, this chapter will focus on the performance of these techniques and tools for plasmonic structures. This is done through an overview of benchmarks available in literature and by considering a few thoroughly analyzed structures.

Missing links will be pointed out and suggestions will be given for the future.

Guy A. E. Vandenbosch

http://dx.doi.org/10.5772/51166

design as real measurements.

**1. Introduction** 

to be solved.

Additional information is available at the end of the chapter


## **Computational Electromagnetics in Plasmonics**

Guy A. E. Vandenbosch

22 Plasmonics – Principles and Applications

2005; 71: 125201-125210.

2165.

9(9): 3349–3353.

[26] Hendry E, Koeberg M, Schins JM, Nienhuys HK, Sundström V, Siebbeles LDA, Bonn M. Interchain effects in the ultrafast photophysics of a semiconducting polymer: THz timedomain spectroscopy of thin films and isolated chains in solution. Physical Review B

[27] Parkinson P, Lloyd-Hughes J, Gao Q, Tan HH, Jagadish C, Johnston MB, Herz LM. Transient Terahertz Conductivity of GaAs Nanowires. Nano Letters 2007; 7(7): 2162–

[28] Parkinson P, Joyce HJ, Gao Q, Tan HH, Zhang X, Zou J, Jagadish C, Herz LM, Johnston MB. Carrier Lifetime and Mobility Enhancement in Nearly Defect-Free Core−Shell Nanowires Measured Using Time-Resolved Terahertz Spectroscopy. Nano Letters 2009;

[29] Di Sia P, Dallacasa V, Dallacasa F. Transient conductivity in nanostructured films**.** 

[30] Di Sia P, Dallacasa V. Anomalous charge transport: a new "time domain"

Journal of Nanoscience and Nanotechnology 2011; 11: 1-6.

generalization of the Drude model. Plasmonics 2011; 6(1): 99-104.

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/51166

### **1. Introduction**

In electromagnetics, numerical techniques have been essential in the development of new technology in the last two decades. The rapidly growing computer capacity and calculation speeds make accurate solutions of very complex problems feasible. This has been especially true in the design of microwave and millimeter wave components and antennas. Whereas 30 years ago, the design of an antenna was based on simple analytical models, or trial and error strategies, nowadays, simulations seem to be as crucial to the design as real measurements.

The situation is quite different in plasmonics. Plasmonics is a quite novel research field and the application of computational electromagnetics in plasmonics can be categorized as "very recent". There are many challenges that still need to be faced and "missing links" that have to be solved.

The plasmonic structures targeted are structures in the order of magnitude of a wavelength at plasmonic frequencies, i.e. at near IR and optical frequencies, and beyond. Although this frequency range is totally different from the traditional range where computational tools have been developed, i.e. the microwave range, in most cases no special modeling techniques have to be used. By far most plasmonic topologies reported in literature have been analyzed / designed with the well-known numerical techniques implemented within in-house developed or commercial software packages. This means that in this chapter the major numerical techniques can be overviewed in a general sense, referring to standard literature. These techniques will not be derived or explained here in full detail. Instead, this chapter focuses on those aspects that come into the picture when the structure is plasmonic.

After the section on techniques and tools, this chapter will focus on the performance of these techniques and tools for plasmonic structures. This is done through an overview of benchmarks available in literature and by considering a few thoroughly analyzed structures. Missing links will be pointed out and suggestions will be given for the future.

© 2012 Vandenbosch, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### **2. Fundamental physical modeling differences**

This section discusses the differences between classical topologies handled with computational electromagnetics and plasmonic topologies. Classical means structures at much lower frequencies, for example at microwave frequencies, where computational electromagnetics is a well-developed mature field.

Computational Electromagnetics in Plasmonics 25

which corresponds to conductivity. These data were obtained through experimental ellipsometry. Note that the difference between the two metals is enormous. It is evident that this strong variation has a serious impact on the behavior of a device over the frequency band. The variability of material characteristics thus has to be taken into account fully into

permittivity of Cu

❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜

200 300 400 500 600 700 800 900 1000 freq. [THz]

permittivity of Ag

❜

❜ ❜

Re(*<sup>m</sup>*) 10 x Im(*<sup>m</sup>*)

the modeling.

**Figure 1.** Top: permittivity of Cu; bottom: permittivity of Ag.



0

20

40

❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜

method in the case of plasmonics is discussed.

In this section the full wave modeling techniques are introduced and categorized on the basis of their solution method: Finite Elements (FE), Finite Differences in the Time Domain (FDTD), Finite Integration Technique (FIT), and Integral Equations (IE) solved by the Method of Moments (MoM). Based on their theoretical specificities, the application of each

200 300 400 500 600 700 800 900 1000 freq. [THz]

❜

❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜ ❜

Re(*<sup>m</sup>*) 10 x Im(*<sup>m</sup>*)

**3. Modeling techniques** 



0

20

❜ ❜

40

### **2.1. Frequency**

First of all, it is essential to point out that the interaction between light and plasmonic structures in the frequency bands considered can still be analyzed with a high degree of accuracy using classical electromagnetic theory. Although the frequency is orders of magnitude higher in plasmonics compared to microwaves, Mawell's laws are the same. They are linear, and thus scalable. As long as quantum effects do not have to be taken into account, which is still the case for the plasmonic applications considered, since the structures are not that small [1], there is no fundamental problem. The underlying formulation of Maxwell's equations can remain unaffected.

The fact that at this small scale, no quantum effects have to be taken into account is really a crucial observation. It means that the concept of a "scatterer", a device able to be excited by electromagnetic waves rather than particles, still works. Basically, the coupling between an electromagnetic (light-) wave and a plasmonic scatterer is thus the same as it is at microwave frequencies, and can be studied in the same way.

### **2.2. Volumetric currents**

A consequence of the high frequencies is the small scale: the elementary building blocks of the topologies considered are at nanoscale. There are two important issues related to this. First, nanoscale fabrication technology of today is only able to generate 3D type structures (i.e. volumes). Thin 2D sheets, where the thickness of the metal is orders of magnitude smaller than the transversal dimensions of the pattern, as commonly used at microwave frequencies (for example a conducting strip or patch), are not possible. Second, for nanostructures operating at plasmonic frequencies the skin depth may be comparable with the structural dimensions, so that currents do flow over the complete volume and the device has to be described with volumetric currents. A surface current description is not sufficient.

### **2.3. Material characteristics**

At microwave frequencies, in most cases the material properties are constant over the frequency bands considered. Also, apart from the losses, most metals behave more or less in the same way, i.e. as good conductors. At IR and optical frequencies however, most metals have very dispersive properties. Permittivities and conductivities may vary orders of magnitude over these frequency bands. The real part of the permittivity may even be negative. In Fig. 1, the permittivity of Cu and Ag are depicted, both real and imaginary part, which corresponds to conductivity. These data were obtained through experimental ellipsometry. Note that the difference between the two metals is enormous. It is evident that this strong variation has a serious impact on the behavior of a device over the frequency band. The variability of material characteristics thus has to be taken into account fully into the modeling.

**Figure 1.** Top: permittivity of Cu; bottom: permittivity of Ag.

#### **3. Modeling techniques**

24 Plasmonics – Principles and Applications

**2.1. Frequency** 

**2.2. Volumetric currents** 

**2.3. Material characteristics** 

**2. Fundamental physical modeling differences** 

formulation of Maxwell's equations can remain unaffected.

microwave frequencies, and can be studied in the same way.

electromagnetics is a well-developed mature field.

This section discusses the differences between classical topologies handled with computational electromagnetics and plasmonic topologies. Classical means structures at much lower frequencies, for example at microwave frequencies, where computational

First of all, it is essential to point out that the interaction between light and plasmonic structures in the frequency bands considered can still be analyzed with a high degree of accuracy using classical electromagnetic theory. Although the frequency is orders of magnitude higher in plasmonics compared to microwaves, Mawell's laws are the same. They are linear, and thus scalable. As long as quantum effects do not have to be taken into account, which is still the case for the plasmonic applications considered, since the structures are not that small [1], there is no fundamental problem. The underlying

The fact that at this small scale, no quantum effects have to be taken into account is really a crucial observation. It means that the concept of a "scatterer", a device able to be excited by electromagnetic waves rather than particles, still works. Basically, the coupling between an electromagnetic (light-) wave and a plasmonic scatterer is thus the same as it is at

A consequence of the high frequencies is the small scale: the elementary building blocks of the topologies considered are at nanoscale. There are two important issues related to this. First, nanoscale fabrication technology of today is only able to generate 3D type structures (i.e. volumes). Thin 2D sheets, where the thickness of the metal is orders of magnitude smaller than the transversal dimensions of the pattern, as commonly used at microwave frequencies (for example a conducting strip or patch), are not possible. Second, for nanostructures operating at plasmonic frequencies the skin depth may be comparable with the structural dimensions, so that currents do flow over the complete volume and the device has to be described with volumetric currents. A surface current description is not sufficient.

At microwave frequencies, in most cases the material properties are constant over the frequency bands considered. Also, apart from the losses, most metals behave more or less in the same way, i.e. as good conductors. At IR and optical frequencies however, most metals have very dispersive properties. Permittivities and conductivities may vary orders of magnitude over these frequency bands. The real part of the permittivity may even be negative. In Fig. 1, the permittivity of Cu and Ag are depicted, both real and imaginary part,

In this section the full wave modeling techniques are introduced and categorized on the basis of their solution method: Finite Elements (FE), Finite Differences in the Time Domain (FDTD), Finite Integration Technique (FIT), and Integral Equations (IE) solved by the Method of Moments (MoM). Based on their theoretical specificities, the application of each method in the case of plasmonics is discussed.

The cradle of computational electromagnetics can be found in the microwave research community. Since this community traditionally is dealing with structures in the order of wavelengths, right from its beginning days, it had no choice than to try to rigorously solve Maxwell's equations. Most standard reference works on full-wave computational electromagnetics by consequence can be found within this community. A history and a comprehensive overview of the different numerical techniques and of their application in computational electromagnetics (CEM) may be found in [2]. Recent work on the last developments in CEM [3] concentrates on the two main approaches of differential and integral methods. A Good review and perspectives concerning the relationship between differential and integral equations (IE) modeling is recommended in a review paper by Miller [4**]**. A thorough discussion on the different techniques with clarifying examples is also given in **[**5].

#### **3.1. Differential equation techniques**

In time domain, Maxwell's equations are

$$
\nabla \times \mathbf{E} = -\frac{d\mathbf{B}}{dt} \tag{1}
$$

with ω π

approximate numerical way.

*3.1.1. The Finite Element Method (FEM) [7], [8]* 

= ⋅ 2 frequency the pulsation. Note that in (5) and (6)

these concepts directly. As will be shown later, this complicates things.

are related to the different ways in which space (and time) can be chopped up.

here they are depending on frequency, and in this way dispersion is fully taken into account (see also section 2.3). In the field of plasmonics, this is an important difference between computational tools in time and in frequency domain. Computational tools in frequency domain may use directly the well-known concepts of permittivity and permeability, albeit that they become frequency dependent. Computational tools in time domain cannot use

From a mathematical perspective, Maxwell's equations are differential equations relating vector fields to each other. Differential equation methods in Computational Electromagnetics are methods that directly consider Maxwell's equations (or the Helmholtz wave equations derived from them), with little analytical preprocessing. Basically, these differential equations, which are valid in any point of space and time, are solved by approximating them by difference equations, which are valid in a discrete set of points in space and time. This is done by chopping up space (and time if time domain is considered) in little pieces in which the field variation has a pre-described profile. This reduces the problem of fields varying over space (and time) to a discrete (matrix) problem that can be handled on a computer. The differences between the several differential equation methods

Since the number of unknowns is proportional to the volume and the resolution considered, differential equation methods are particularly suitable for modeling small full threedimensional volumes that have complex geometrical details, for example smaller closedregion problems involving inhomogeneous media [6]. Intrinsically, differential equations are less suited for open problems. The reason is that in principle they require a discretization of the entire space under consideration. This space is limited in case of closed problems, but corresponds to infinite space in case of open problems. In practice, this problem is solved by the introduction of techniques like Absorbing Boundary Conditions, and Perfectly Matched Layers (PML) [7]. They mimic the wave arriving at the boundary as travelling further to infinity. The quality of these truncating techniques nowadays is very high so that, in practice, the intrinsic problem with open structures has been overcome, albeit it in an

The most popular differential equation-based methods are the Finite Element Method (FEM), for example utilized in Ansoft's HFSS software package, and the Finite-Difference Time Domain method (FDTD), which is employed for example by CST's Time Domain

FEM is a method based on solving partial differential equations. It is most commonly formulated based on a variational expression. It subdivides space in elements, for example tetrahedra. Fields inside these elements are expressed in terms of a number of basic functions, for example polynomials. These expressions are inserted into the functional of the equations, and the variation of the functional is made zero. This yields a matrix eigenvalue

transient solver (in the particular case of Cartesian grids), and by Lumerical.

Computational Electromagnetics in Plasmonics 27

are used. In general,

ε and μ

$$\nabla \times \mathbf{H} = \frac{d\mathbf{D}}{dt} + \mathbf{J} \tag{2}$$

with **E** and **B** the electric field and the magnetic induction, respectively, **H** and **D** the magnetic field and the electric induction, respectively, and **J** the electric current flowing. In free space, and in a homogeneous, isotropic, time-invariant, linear medium (most materials behave like this in the microwave frequency range), the following relations hold

$$\mathbf{D} = \mathbf{\mathcal{E}} \tag{3}$$

$$\mathbf{B} = \mu \mathbf{H} \tag{4}$$

They are called the constitutive relations, with ε and μ the permittivity and permeability of the medium surrounding the observation point considered. It is crucial to point out that in plasmonics in general (3) and (4) cannot be used. For example, the very dispersive properties of metals at optical frequencies (and beyond) prohibit this for these materials.

In time domain, the time variation has to be determined. In frequency domain, it is assumed that all field quantities are varying in a sinusoidal way. This means that the variation with time is known, and only the variation of the fields in space has to be determined. Using complex notation for the frequency domain with a exp( ) *j t* ω time dependency, Maxwell's equations become

$$
\nabla \times \mathbf{E} = -j\alpha \mathbf{B} \tag{5}
$$

$$\nabla \times (\mu^{-1} \mathbf{B}) = j\alpha \mathbf{e} \mathbf{E} + \mathbf{J} \tag{6}$$

with ω π = ⋅ 2 frequency the pulsation. Note that in (5) and (6) ε and μ are used. In general, here they are depending on frequency, and in this way dispersion is fully taken into account (see also section 2.3). In the field of plasmonics, this is an important difference between computational tools in time and in frequency domain. Computational tools in frequency domain may use directly the well-known concepts of permittivity and permeability, albeit that they become frequency dependent. Computational tools in time domain cannot use these concepts directly. As will be shown later, this complicates things.

From a mathematical perspective, Maxwell's equations are differential equations relating vector fields to each other. Differential equation methods in Computational Electromagnetics are methods that directly consider Maxwell's equations (or the Helmholtz wave equations derived from them), with little analytical preprocessing. Basically, these differential equations, which are valid in any point of space and time, are solved by approximating them by difference equations, which are valid in a discrete set of points in space and time. This is done by chopping up space (and time if time domain is considered) in little pieces in which the field variation has a pre-described profile. This reduces the problem of fields varying over space (and time) to a discrete (matrix) problem that can be handled on a computer. The differences between the several differential equation methods are related to the different ways in which space (and time) can be chopped up.

Since the number of unknowns is proportional to the volume and the resolution considered, differential equation methods are particularly suitable for modeling small full threedimensional volumes that have complex geometrical details, for example smaller closedregion problems involving inhomogeneous media [6]. Intrinsically, differential equations are less suited for open problems. The reason is that in principle they require a discretization of the entire space under consideration. This space is limited in case of closed problems, but corresponds to infinite space in case of open problems. In practice, this problem is solved by the introduction of techniques like Absorbing Boundary Conditions, and Perfectly Matched Layers (PML) [7]. They mimic the wave arriving at the boundary as travelling further to infinity. The quality of these truncating techniques nowadays is very high so that, in practice, the intrinsic problem with open structures has been overcome, albeit it in an approximate numerical way.

The most popular differential equation-based methods are the Finite Element Method (FEM), for example utilized in Ansoft's HFSS software package, and the Finite-Difference Time Domain method (FDTD), which is employed for example by CST's Time Domain transient solver (in the particular case of Cartesian grids), and by Lumerical.

### *3.1.1. The Finite Element Method (FEM) [7], [8]*

26 Plasmonics – Principles and Applications

**3.1. Differential equation techniques** 

In time domain, Maxwell's equations are

They are called the constitutive relations, with

equations become

complex notation for the frequency domain with a exp( ) *j t*

given in **[**5].

The cradle of computational electromagnetics can be found in the microwave research community. Since this community traditionally is dealing with structures in the order of wavelengths, right from its beginning days, it had no choice than to try to rigorously solve Maxwell's equations. Most standard reference works on full-wave computational electromagnetics by consequence can be found within this community. A history and a comprehensive overview of the different numerical techniques and of their application in computational electromagnetics (CEM) may be found in [2]. Recent work on the last developments in CEM [3] concentrates on the two main approaches of differential and integral methods. A Good review and perspectives concerning the relationship between differential and integral equations (IE) modeling is recommended in a review paper by Miller [4**]**. A thorough discussion on the different techniques with clarifying examples is also

> *d dt*

*d dt*

with **E** and **B** the electric field and the magnetic induction, respectively, **H** and **D** the magnetic field and the electric induction, respectively, and **J** the electric current flowing. In free space, and in a homogeneous, isotropic, time-invariant, linear medium (most materials

> **D E** = ε

**B H** = μ

of the medium surrounding the observation point considered. It is crucial to point out that in plasmonics in general (3) and (4) cannot be used. For example, the very dispersive properties of metals at optical frequencies (and beyond) prohibit this for these materials.

In time domain, the time variation has to be determined. In frequency domain, it is assumed that all field quantities are varying in a sinusoidal way. This means that the variation with time is known, and only the variation of the fields in space has to be determined. Using

> ∇× =− **E B** *j*ω

> > ωε

<sup>1</sup> ( ) μ

ω

*j* <sup>−</sup> ∇× = + **B EJ** (6)

ε and μ

behave like this in the microwave frequency range), the following relations hold

∇× =− **<sup>B</sup> <sup>E</sup>** (1)

∇× = + **<sup>D</sup> H J** (2)

(3)

(4)

the permittivity and permeability

time dependency, Maxwell's

(5)

FEM is a method based on solving partial differential equations. It is most commonly formulated based on a variational expression. It subdivides space in elements, for example tetrahedra. Fields inside these elements are expressed in terms of a number of basic functions, for example polynomials. These expressions are inserted into the functional of the equations, and the variation of the functional is made zero. This yields a matrix eigenvalue

#### 28 Plasmonics – Principles and Applications

equation whose solution yields the fields at the nodes. FEM gives rise to a very sparse matrix equation, which can be solved using dedicated matrix algebra technology, leading to very fast solution times, considering the huge number of unknowns.

Computational Electromagnetics in Plasmonics 29

FDTD technique dedicated to plasmonics and photonics can be found in the commercial tool

The Finite Integration Technique was introduced by Weiland in 1977. The word integration does not imply any relation with integral equations. FIT first describes Maxwell's equations on a grid space. The matrix equations for the electromagnetic integral quantities obtained by FIT possess some inherent properties of Maxwell's equations, for example with respect to charge and energy conservation. This makes them very attractive from a theoretical point of view. FIT can be formulated on different kinds of grids, e.g. Cartesian or general nonorthogonal ones, which is a clear advantage. In the time-domain, the resulting discrete grid equations of FIT are, at least in "some" cases, identical to the discrete equations derived with the classical Finite-Difference Time-Domain (FDTD) method. In contrast to FIT, which is applied to the integral form of the field equations, FDTD (as a subset of the finite integration method) is applied to the differential form of the governing Maxwell curl equations. Theoretical links between FDTD and the FIT approach realized in CST may be found in [17]. A comparative study on Time Domain methods was presented in [18]. In some sense, FIT can thus be considered as a powerful generalization of the FDTD technique. A software tool

Integral equation methods make use of Maxwell's equations in integral equation form to formulate the electromagnetic problem in terms of unknown currents flowing on the object to be described. These currents are induced by a field incident on the object. This incident field can be a real incident field, travelling in space, or a bounded wave feeding the object for example via a transmission line. An integral equation solution is fundamentally based on

**Figure 2.** Standard Cartesian Yee cell used in the FDTD technique.

using FIT and very widely spread is CST Microwave Studio.

**3.2. Integral Equation techniques (IE)** 

the combination of two equations.

*3.1.3. The Finite Integration Technique (FIT) [16]* 

Lumerical [15].

Its first formulations were developed as matrix methods for structural mechanics. This lead to the idea to approximate solids and Courant (1942) introduced an assembly of triangular elements and the minimum of potential energy to torsion problems [9]. The first paper on the application of FEM to electrical problems appeared in 1968 [10]. An extensive review on the history of FEM in electromagnetics was published in an issue of the Antennas and Propagation Magazine [11]. FEM normally is formulated in the frequency domain, i.e. for time-harmonic problems. This means that, as for IE-MoM, the solution has to be calculated for every frequency of interest.

Numerous references can be given developing, explaining, and using FEM. A good book to start with is [8]. A software tool using FEM and very widely spread is Ansoft HFSS.

### *3.1.2. The Finite-Difference Time-Domain technique (FDTD) [12], [13], [14]*

The nature of Maxwell's differential equations is that the time derivative of the H-field is dependent on the curl of the E-field, and the time derivative of the H-field is dependent on the curl of the E-field. These basic properties result in the core FDTD time-stepping relation that, at any point in space, an updated value of an E/H-field in time is dependent on the stored value of the E/H-field and the numerical curl of the local distribution of the H/E-field in space. The numerical translation into a time-stepping algorithm was introduced by Yee in 1966. Indeed, swapping between E-field and H-field updates allows to define a marchingon-in-time process wherein sampled fields of the continuous electromagnetic waves under consideration are used. These waves can be seen to propagate in the Yee lattice, a numerical three-dimensional space lattice comprised of a multiplicity of Yee cells, see Fig. 2. More specifically, Yee proposed a leapfrog scheme for marching-on in time wherein the E-field and H-field updates are staggered so that E-field updates are observed midway during each time-step between successive H-field updates, and vice versa. A huge advantage of FDTD is that this explicit time-stepping scheme avoids the need to solve simultaneous equations, so no matrix inversions are necessary. It also yields dissipation-free numerical wave propagation. Negative is that this scheme results in an upper bound on the time-step to ensure numerical stability. This means that simulations may require many thousands of time-steps for completion. The use of the Yee lattice has proven to be very robust in numerical calculations.

FDTD is extremely versatile since the interaction of an electromagnetic wave with matter can be mapped into the space lattice by assigning appropriate values of permittivity to each electric field component, and permeability to each magnetic field component. This can be done without seriously compromising the speed of the method.

The fact that time is observed directly, and the versatility of the method, make it probably the most efficient technique for complex 3D transient problems. An implementation of the FDTD technique dedicated to plasmonics and photonics can be found in the commercial tool Lumerical [15].

**Figure 2.** Standard Cartesian Yee cell used in the FDTD technique.

#### *3.1.3. The Finite Integration Technique (FIT) [16]*

28 Plasmonics – Principles and Applications

for every frequency of interest.

numerical calculations.

equation whose solution yields the fields at the nodes. FEM gives rise to a very sparse matrix equation, which can be solved using dedicated matrix algebra technology, leading to

Its first formulations were developed as matrix methods for structural mechanics. This lead to the idea to approximate solids and Courant (1942) introduced an assembly of triangular elements and the minimum of potential energy to torsion problems [9]. The first paper on the application of FEM to electrical problems appeared in 1968 [10]. An extensive review on the history of FEM in electromagnetics was published in an issue of the Antennas and Propagation Magazine [11]. FEM normally is formulated in the frequency domain, i.e. for time-harmonic problems. This means that, as for IE-MoM, the solution has to be calculated

Numerous references can be given developing, explaining, and using FEM. A good book to

The nature of Maxwell's differential equations is that the time derivative of the H-field is dependent on the curl of the E-field, and the time derivative of the H-field is dependent on the curl of the E-field. These basic properties result in the core FDTD time-stepping relation that, at any point in space, an updated value of an E/H-field in time is dependent on the stored value of the E/H-field and the numerical curl of the local distribution of the H/E-field in space. The numerical translation into a time-stepping algorithm was introduced by Yee in 1966. Indeed, swapping between E-field and H-field updates allows to define a marchingon-in-time process wherein sampled fields of the continuous electromagnetic waves under consideration are used. These waves can be seen to propagate in the Yee lattice, a numerical three-dimensional space lattice comprised of a multiplicity of Yee cells, see Fig. 2. More specifically, Yee proposed a leapfrog scheme for marching-on in time wherein the E-field and H-field updates are staggered so that E-field updates are observed midway during each time-step between successive H-field updates, and vice versa. A huge advantage of FDTD is that this explicit time-stepping scheme avoids the need to solve simultaneous equations, so no matrix inversions are necessary. It also yields dissipation-free numerical wave propagation. Negative is that this scheme results in an upper bound on the time-step to ensure numerical stability. This means that simulations may require many thousands of time-steps for completion. The use of the Yee lattice has proven to be very robust in

FDTD is extremely versatile since the interaction of an electromagnetic wave with matter can be mapped into the space lattice by assigning appropriate values of permittivity to each electric field component, and permeability to each magnetic field component. This can be

The fact that time is observed directly, and the versatility of the method, make it probably the most efficient technique for complex 3D transient problems. An implementation of the

done without seriously compromising the speed of the method.

start with is [8]. A software tool using FEM and very widely spread is Ansoft HFSS.

*3.1.2. The Finite-Difference Time-Domain technique (FDTD) [12], [13], [14]* 

very fast solution times, considering the huge number of unknowns.

The Finite Integration Technique was introduced by Weiland in 1977. The word integration does not imply any relation with integral equations. FIT first describes Maxwell's equations on a grid space. The matrix equations for the electromagnetic integral quantities obtained by FIT possess some inherent properties of Maxwell's equations, for example with respect to charge and energy conservation. This makes them very attractive from a theoretical point of view. FIT can be formulated on different kinds of grids, e.g. Cartesian or general nonorthogonal ones, which is a clear advantage. In the time-domain, the resulting discrete grid equations of FIT are, at least in "some" cases, identical to the discrete equations derived with the classical Finite-Difference Time-Domain (FDTD) method. In contrast to FIT, which is applied to the integral form of the field equations, FDTD (as a subset of the finite integration method) is applied to the differential form of the governing Maxwell curl equations. Theoretical links between FDTD and the FIT approach realized in CST may be found in [17]. A comparative study on Time Domain methods was presented in [18]. In some sense, FIT can thus be considered as a powerful generalization of the FDTD technique. A software tool using FIT and very widely spread is CST Microwave Studio.

#### **3.2. Integral Equation techniques (IE)**

Integral equation methods make use of Maxwell's equations in integral equation form to formulate the electromagnetic problem in terms of unknown currents flowing on the object to be described. These currents are induced by a field incident on the object. This incident field can be a real incident field, travelling in space, or a bounded wave feeding the object for example via a transmission line. An integral equation solution is fundamentally based on the combination of two equations.

In the case of Electric Field Integral Equations (the case of prime importance to the plasmonics community), the solution is based on considering the electric field. The first equation is

$$\mathbf{E}^{\text{sca}}(\mathbf{r}) = \int\_{V'} \mathbf{G}(\mathbf{r}, \mathbf{r}') \cdot \mathbf{J}^{ind}(\mathbf{r}') dV' \tag{7}$$

Computational Electromagnetics in Plasmonics 31

Next to the reduced number of unknowns, there are other theoretical advantages linked to the IE-MoM method. The second advantage is that, if properly formulated, IE-MoM *is variationally stable* since most of the output parameters are expressed in integral form over the equivalent currents. This means that even if the calculated currents differ considerably from the exact solution, integral parameters over both currents may remain very similar. Further, this will be illustrated by showing that even with a rather rough mesh high quality physical results may be obtained. A third advantage is that MoM does not heavily suffer from field singularities for example near sharp edges, since they are analytically incorporated inside the Green's functions. For differential equation methods special care (= a

The main disadvantage of IE-MoM is that, although there are many efforts in that direction, there is still a lack of matrix solvers operating on dense matrices which are comparable in efficiency to the solvers used in the differential equation techniques (mainly FEM), which yield sparse matrices. This has precluded the use of the very flexible volumetric IE-MoM

IE-MoM is normally applied in the frequency domain, i.e. for time-harmonic problems. This

Much more details can be found in the classic and basic Method of Moments (MoM) book [19]. There are many variants of the method. In general, boundary equations can be enforced also at the boundaries of volumes (utilizing Surface Integral Equations (SIE) [20]), next to inside the entire volumes themselves (applying Volume Integral Equations (VIE) [21] at the inside of the components, as described above). Also, the integral equations can be written down in different forms (dyadic form, mixed-potential form, hybrid forms, etc. [22]), which give rise to specific implementations. Further, the first theoretical developments in the field

In the following sections, several aspects linked with the use of integral equation techniques

The fact that a 3D volumetric current has to be described, as mentioned already, invokes the need for either a real volumetric MoM implementation or the more common surface approach, where the plasmonic component is described by applying the equivalence principle at its surface. The surface approach has already been used, for example in [24]. The volumetric approach as described in [21], has very recently been introduced in the plasmonics community [25]. For realistic structures, where volumes can be embedded

Although strongly varying material characteristics may seem trivial to implement in an integral equation scheme, in some cases more severe consequences occur. For example,

fine mesh) should be taken to describe correctly these field singularities.

technique for modeling complex structures on widespread computer systems.

means that the solution has to be determined at each frequency of interest.

of computational plasmonics are appearing in literature [23].

within a layer structure, it offers a very flexible solution technique.

for plasmonic structures are discussed.

*3.2.1. 3D volumetric current* 

*3.2.2. Material characteristics* 

which gives the scattered electric field *sca* **E** generated in an arbitrary observation point **r** in space in terms of an induced volumetric electric current distribution *ind* **J** flowing within the volume *V* ' of the object considered. The kernel **Grr** ( , ') is a so-called dyadic Green's function, which is a tensor, and which relates a Dirac impulse type current flowing in **r**' to the electric field it generates. The main strength of the integral equation technique is that in many cases the dyadic Green's function can be calculated either directly analytically, such as in homogeneous media, or as some kind of inverse Fourier transform of an analytic function, as is the case for example in layered media. It is crucial to understand that the dyadic Green's function is actually a solution of Maxwell's equations and thus rigorously takes into account the background medium, for example a stack of dielectric layers. This means that unknown currents only have to be assumed on the objects embedded within this background medium. This results in an enormous reduction of unknowns compared to differential equation techniques, which have to model these dielectric layers in just the same way as the objects embedded within them. Basically, the problem formulation automatically covers the entire surrounding space without making any fundamental approximations. As a consequence, the corresponding solution is automatically valid in every point of the background medium. Far field radiation phenomena, surface waves in layered structures, etc., that are vital for efficient and accurate analysis, are analytically included in the solution.

In the case of plasmonics, the second equation is

$$\begin{aligned} \mathbf{J}^{ind}(\mathbf{r}') &= j\phi \left( \varepsilon(\mathbf{r}') - \varepsilon\_0 \right) \mathbf{E}^{ind}(\mathbf{r}')\\ &= j\phi \left( \varepsilon(\mathbf{r}') - \varepsilon\_0 \right) \left( \mathbf{E}^{sca}(\mathbf{r}') + \mathbf{E}^{sca}(\mathbf{r}') \right) \end{aligned} \tag{8}$$

It expresses the boundary condition which has to be enforced within the object under consideration. It is a relation between the total field, which is the sum of the incident and the scattered field, and the volumetric electric current, which can be considered as partially being a conduction current (due to the imaginary part of the permittivity), and a polarization current (due to the real part of the permittivity).

Combining (7) and (8) yields the equation from which the currents can be solved. This can be done by chopping up the object considered (only the object, and thus not entire space !!!) in little pieces in which the current variation has a pre-described profile (a so-called basis function). This reduces the problem of currents varying over the object to a discrete (matrix) problem that can be handled on a computer. IE-MoM gives rise to a dense matrix equation, which can be solved using standard matrix algebra technology.

Next to the reduced number of unknowns, there are other theoretical advantages linked to the IE-MoM method. The second advantage is that, if properly formulated, IE-MoM *is variationally stable* since most of the output parameters are expressed in integral form over the equivalent currents. This means that even if the calculated currents differ considerably from the exact solution, integral parameters over both currents may remain very similar. Further, this will be illustrated by showing that even with a rather rough mesh high quality physical results may be obtained. A third advantage is that MoM does not heavily suffer from field singularities for example near sharp edges, since they are analytically incorporated inside the Green's functions. For differential equation methods special care (= a fine mesh) should be taken to describe correctly these field singularities.

The main disadvantage of IE-MoM is that, although there are many efforts in that direction, there is still a lack of matrix solvers operating on dense matrices which are comparable in efficiency to the solvers used in the differential equation techniques (mainly FEM), which yield sparse matrices. This has precluded the use of the very flexible volumetric IE-MoM technique for modeling complex structures on widespread computer systems.

IE-MoM is normally applied in the frequency domain, i.e. for time-harmonic problems. This means that the solution has to be determined at each frequency of interest.

Much more details can be found in the classic and basic Method of Moments (MoM) book [19]. There are many variants of the method. In general, boundary equations can be enforced also at the boundaries of volumes (utilizing Surface Integral Equations (SIE) [20]), next to inside the entire volumes themselves (applying Volume Integral Equations (VIE) [21] at the inside of the components, as described above). Also, the integral equations can be written down in different forms (dyadic form, mixed-potential form, hybrid forms, etc. [22]), which give rise to specific implementations. Further, the first theoretical developments in the field of computational plasmonics are appearing in literature [23].

In the following sections, several aspects linked with the use of integral equation techniques for plasmonic structures are discussed.

### *3.2.1. 3D volumetric current*

30 Plasmonics – Principles and Applications

equation is

In the case of Electric Field Integral Equations (the case of prime importance to the plasmonics community), the solution is based on considering the electric field. The first

( ) ( , ') ( ') ' *sca ind*

which gives the scattered electric field *sca* **E** generated in an arbitrary observation point **r** in

volume *V* ' of the object considered. The kernel **Grr** ( , ') is a so-called dyadic Green's function, which is a tensor, and which relates a Dirac impulse type current flowing in **r**' to the electric field it generates. The main strength of the integral equation technique is that in many cases the dyadic Green's function can be calculated either directly analytically, such as in homogeneous media, or as some kind of inverse Fourier transform of an analytic function, as is the case for example in layered media. It is crucial to understand that the dyadic Green's function is actually a solution of Maxwell's equations and thus rigorously takes into account the background medium, for example a stack of dielectric layers. This means that unknown currents only have to be assumed on the objects embedded within this background medium. This results in an enormous reduction of unknowns compared to differential equation techniques, which have to model these dielectric layers in just the same way as the objects embedded within them. Basically, the problem formulation automatically covers the entire surrounding space without making any fundamental approximations. As a consequence, the corresponding solution is automatically valid in every point of the background medium. Far field radiation phenomena, surface waves in layered structures, etc., that are vital for efficient and accurate analysis, are analytically included in the solution.

> 0 0

 ε

 ε

=− +

*sca sca*

**r Er Er** (8)

( ( ') ) ( ( ') ( '))

It expresses the boundary condition which has to be enforced within the object under consideration. It is a relation between the total field, which is the sum of the incident and the scattered field, and the volumetric electric current, which can be considered as partially being a conduction current (due to the imaginary part of the permittivity), and a

Combining (7) and (8) yields the equation from which the currents can be solved. This can be done by chopping up the object considered (only the object, and thus not entire space !!!) in little pieces in which the current variation has a pre-described profile (a so-called basis function). This reduces the problem of currents varying over the object to a discrete (matrix) problem that can be handled on a computer. IE-MoM gives rise to a dense matrix equation,

( ') ( ( ') ) ( ')

*ind ind*

**Jr r Er**

= −

ωε

*j j* ωε

polarization current (due to the real part of the permittivity).

which can be solved using standard matrix algebra technology.

= ⋅ *dV* **E r Grr J r** (7)

**J** flowing within the

'

*V*

space in terms of an induced volumetric electric current distribution *ind*

In the case of plasmonics, the second equation is

The fact that a 3D volumetric current has to be described, as mentioned already, invokes the need for either a real volumetric MoM implementation or the more common surface approach, where the plasmonic component is described by applying the equivalence principle at its surface. The surface approach has already been used, for example in [24]. The volumetric approach as described in [21], has very recently been introduced in the plasmonics community [25]. For realistic structures, where volumes can be embedded within a layer structure, it offers a very flexible solution technique.

#### *3.2.2. Material characteristics*

Although strongly varying material characteristics may seem trivial to implement in an integral equation scheme, in some cases more severe consequences occur. For example,

#### 32 Plasmonics – Principles and Applications

some Green's function calculation schemes are based on the fact that the material properties do not change over the frequency band of interest. Any particular implementation of IE-MoM needs this essential relaxation into its formalism, which possibly needs to be adapted.

Computational Electromagnetics in Plasmonics 33

FEM solution of the electromagnetic topology under consideration. This software is extremely popular and is used for all kinds of purposes. Numerous results for plasmonic

The first step in a new simulation project is to define the geometry of the structure. The geometry is modeled in the GUI or can be imported from another program (AutoCAD, STEP, …). Then the user continues by defining material properties, boundary conditions, and excitations to the different domains and surfaces of the geometry. If desired one can control the meshing process manually. Finally, the solution parameters are defined. The most important are the solution frequency, the order of the base functions and the matrix solver type. The direct matrix solver is default, but an iterative solver can be used as well. Very useful is the automatic adaptive mesh generation and refinement, which in many cases frees the designer of worrying about which mesh/grid to choose. First a solution for an initial mesh is determined and the solution is assessed. If the solution does not qualify, the mesh is refined and a new solution is computed. This procedure is repeated until one of the exit criteria is fulfilled. It is important to note that HFSS offers no curved elements for a better approximation of curved objects. This means that the mesh along a curved surface has to be chosen very fine. This can be achieved with the 'surface approximation' parameter. Finally the solution, i.e. electric and magnetic fields, currents, and S-parameters can be visualized in 1D, 2D and 3D. All solutions can be exported as files. An 'Optimetrics' toolbox

topologies can be easily found by googling the words HFSS and plasmonics.

is offered for optimizations, parameter, sensitivity and statistical analysis.

Comsol Multiphysics uses the finite element method to solve partial differential equations in 2D and 3D. It is a multiphysics code, which means that it can handle not only electromagnetics, but also acoustics, mechanics, fluid dynamics, heat transfer, etc.. The geometry, material parameters and boundary conditions have to be set up in the graphical user interface (GUI). Then, after defining boundaries and domains the mesh is generated. This includes the selection of the order of the basis functions and the curved mesh elements (order higher than one means curved element). They are used for a better approximation of curved boundaries. The meshing can be steered by specifying a number of nodes and their distribution on each edge of the model. Comsol Multiphysics offers a wide variety of matrix solvers to solve the system of equations. There are direct solvers, more suitable for smaller problems, and iterative solvers suitable for larger problems. Comsol Multiphysics can generate 2D and 3D plots. Quantities derived from the electromagnetic field, like Poynting vector, energy, and energy loss, are available as predefined variables. Comsol Multiphysics does not allow parameter sweeps or optimization. COMSOL is very easily combined with MATLAB. A simulation set up with the GUI can be saved as a Matlab script, which can be easily modified. The Matlab optimization toolbox can be used to optimize a topology.

JCMsuite is a software package based on FEM. It contains dedicated tools for certain simulation problems appearing in nano-optics and plasmonics. It incorporates scattering

*4.1.2. COMSOL multiphysics [28]: FEM* 

*4.1.3. JCMsuite [29]: FEM* 

An advantage of a frequency domain technique like IE-MoM compared to a time domain technique like FDTD is that the ellipsometric measurement data can be used directly in the tool, without any further fitting. FDTD techniques developed for the optical range tend to fit the dielectric response of the metal to the experimentally determined dielectric permittivity using a Drude model [26] or more sophisticated models (for example in Lumerical). However, this may create problems (see further).

### *3.2.3. Occurrence of layer structures*

As explained above, IE technique solvers are formulated making use of Green's functions. These Green's functions can be formulated for multi-layered structures, such that the background medium of the structure may consist of an arbitrary number of horizontal, infinitely stretched, dielectric and metallic layers, which are taken into account analytically. This environment is particularly interesting for plasmonic structures, because they are traditionally deposited on a flat layered dielectric substrate, for example glass. This glass can be contained in the background environment. The only remaining components are local scattering objects with medium or small dimensions compared to the wavelength, such as the dipole in section 5.2.1, but also dots, rods, monomers, dimmers, rings, discs, etc.. With volumetric integral equations these components are replaced by equivalent volume currents, which appear as the primary unknowns in the resulting integral equations. Since the substrates used in normal circumstances are huge compared to the wavelength, edge effects are extremely small and can be neglected. It may thus be concluded that plasmonic structures have specific features which make a IE-MoM solution quite attractive.

### **4. Software tools**

In the following sections, several solvers are briefly described. Several commercial solvers and one academic solver are considered. It has to be emphasized that the author does not claim that this overview is complete. Since the widespread use of computational electromagnetics in plasmonics is quite recent, it is highly probably that there are more solvers that the author is not aware of.

### **4.1. Commercial software tools**

#### *4.1.1. HFSS [27]: FEM*

Since it was one of the first tools in the market, and also due to its generality and flexibility, HFFS is one of the tools heavily used in industrial microwave and millimeter wave design environments. The purpose of HFSS is to extract parasitic parameters (S, Y, Z), visualize 3D electromagnetic fields (near- and far-field), and generate SPICE models, all based on a 3D FEM solution of the electromagnetic topology under consideration. This software is extremely popular and is used for all kinds of purposes. Numerous results for plasmonic topologies can be easily found by googling the words HFSS and plasmonics.

The first step in a new simulation project is to define the geometry of the structure. The geometry is modeled in the GUI or can be imported from another program (AutoCAD, STEP, …). Then the user continues by defining material properties, boundary conditions, and excitations to the different domains and surfaces of the geometry. If desired one can control the meshing process manually. Finally, the solution parameters are defined. The most important are the solution frequency, the order of the base functions and the matrix solver type. The direct matrix solver is default, but an iterative solver can be used as well. Very useful is the automatic adaptive mesh generation and refinement, which in many cases frees the designer of worrying about which mesh/grid to choose. First a solution for an initial mesh is determined and the solution is assessed. If the solution does not qualify, the mesh is refined and a new solution is computed. This procedure is repeated until one of the exit criteria is fulfilled. It is important to note that HFSS offers no curved elements for a better approximation of curved objects. This means that the mesh along a curved surface has to be chosen very fine. This can be achieved with the 'surface approximation' parameter. Finally the solution, i.e. electric and magnetic fields, currents, and S-parameters can be visualized in 1D, 2D and 3D. All solutions can be exported as files. An 'Optimetrics' toolbox is offered for optimizations, parameter, sensitivity and statistical analysis.

### *4.1.2. COMSOL multiphysics [28]: FEM*

32 Plasmonics – Principles and Applications

However, this may create problems (see further).

*3.2.3. Occurrence of layer structures* 

**4. Software tools** 

solvers that the author is not aware of.

**4.1. Commercial software tools** 

*4.1.1. HFSS [27]: FEM* 

some Green's function calculation schemes are based on the fact that the material properties do not change over the frequency band of interest. Any particular implementation of IE-MoM needs this essential relaxation into its formalism, which possibly needs to be adapted. An advantage of a frequency domain technique like IE-MoM compared to a time domain technique like FDTD is that the ellipsometric measurement data can be used directly in the tool, without any further fitting. FDTD techniques developed for the optical range tend to fit the dielectric response of the metal to the experimentally determined dielectric permittivity using a Drude model [26] or more sophisticated models (for example in Lumerical).

As explained above, IE technique solvers are formulated making use of Green's functions. These Green's functions can be formulated for multi-layered structures, such that the background medium of the structure may consist of an arbitrary number of horizontal, infinitely stretched, dielectric and metallic layers, which are taken into account analytically. This environment is particularly interesting for plasmonic structures, because they are traditionally deposited on a flat layered dielectric substrate, for example glass. This glass can be contained in the background environment. The only remaining components are local scattering objects with medium or small dimensions compared to the wavelength, such as the dipole in section 5.2.1, but also dots, rods, monomers, dimmers, rings, discs, etc.. With volumetric integral equations these components are replaced by equivalent volume currents, which appear as the primary unknowns in the resulting integral equations. Since the substrates used in normal circumstances are huge compared to the wavelength, edge effects are extremely small and can be neglected. It may thus be concluded that plasmonic

structures have specific features which make a IE-MoM solution quite attractive.

In the following sections, several solvers are briefly described. Several commercial solvers and one academic solver are considered. It has to be emphasized that the author does not claim that this overview is complete. Since the widespread use of computational electromagnetics in plasmonics is quite recent, it is highly probably that there are more

Since it was one of the first tools in the market, and also due to its generality and flexibility, HFFS is one of the tools heavily used in industrial microwave and millimeter wave design environments. The purpose of HFSS is to extract parasitic parameters (S, Y, Z), visualize 3D electromagnetic fields (near- and far-field), and generate SPICE models, all based on a 3D Comsol Multiphysics uses the finite element method to solve partial differential equations in 2D and 3D. It is a multiphysics code, which means that it can handle not only electromagnetics, but also acoustics, mechanics, fluid dynamics, heat transfer, etc.. The geometry, material parameters and boundary conditions have to be set up in the graphical user interface (GUI). Then, after defining boundaries and domains the mesh is generated. This includes the selection of the order of the basis functions and the curved mesh elements (order higher than one means curved element). They are used for a better approximation of curved boundaries. The meshing can be steered by specifying a number of nodes and their distribution on each edge of the model. Comsol Multiphysics offers a wide variety of matrix solvers to solve the system of equations. There are direct solvers, more suitable for smaller problems, and iterative solvers suitable for larger problems. Comsol Multiphysics can generate 2D and 3D plots. Quantities derived from the electromagnetic field, like Poynting vector, energy, and energy loss, are available as predefined variables. Comsol Multiphysics does not allow parameter sweeps or optimization. COMSOL is very easily combined with MATLAB. A simulation set up with the GUI can be saved as a Matlab script, which can be easily modified. The Matlab optimization toolbox can be used to optimize a topology.

#### *4.1.3. JCMsuite [29]: FEM*

JCMsuite is a software package based on FEM. It contains dedicated tools for certain simulation problems appearing in nano-optics and plasmonics. It incorporates scattering

#### 34 Plasmonics – Principles and Applications

tools, propagation mode tools, and resonance mode tools. The time-harmonic problems can be formulated in 1D, 2D and 3D. The topologies can be isolated or may occur in periodic patterns, or a mixture of both. Interesting is that dedicated tools for problems posed on cylindrically symmetric geometries are available. The electromagnetic fields are discretized with higher order edge elements. JCMsuite contains an automatic mesh generator (with adaptive features), goal-oriented error estimators for adaptive grid refinement, domaindecomposition techniques and fast solvers.

Computational Electromagnetics in Plasmonics 35

time domain solver, it has a sophisticated library of Lorentz, Drude, Debye models for the material parameters. Excitation of the structure is possible with waveguide sources, dipoles, plane waves, focused beams and diffraction-limited spots. A scripting language is available

MAGMAS 3D is the IE-MoM code developed at the Katholieke Universiteit Leuven, Belgium. It was originally developed in cooperation with the European Space Agency for planar and quasi-planar antenna and scattering structures embedded in a multilayered dielectric background medium and operating in the microwave frequency range. Starting from the topology considered, frequency band needed, and type of excitation, it calculates the network, radiation, and scattering characteristics of the structure under consideration. It is based on a full-wave Mixed-Potential formulation of Electric Field Integral Equations (MP-EFIE), originally applied only to 2D surface currents [22], [32]. New theoretical techniques were developed and implemented within the framework: the Expansion Wave Concept (EWC) [33], [34], the Dipole Modeling Technique (DMT) [35], special deembedding procedures [36], etc.. In 2007, it was extended with the capability to handle volumetric 3D currents with the VIE technique introduced in [21]. This was crucial in view of its application in plasmonics. In 2009, this led to the first verified simulated results in this field [25]. Now, it is extensively used for plasmonics. The MAGMAS mesh is based on a combination of rectangular and triangular mesh cells. Full mesh control is available in

To the best knowledge of the author, to date (June 2012), it is the only IE-MoM framework able to handle arbitrary plasmonic structures embedded in a multilayered environment, with the very general and flexible VIE technique. No similar commercial solvers exist, in contrast to the situation at microwave frequencies, where a multitude of IE solvers can be

At this moment plasmonics is a mainly experimentally driven research field. Whereas for example in the field of traditional antenna research, numerous papers can be found on modeling as such, this is much less the case in plasmonics, especially in the high impact journals. Most papers are concerned with the description of the physical phenomena

Most researchers in the plasmonics field use the commercial tools available, and as far as I can see, the main one is Lumerical. In the publications that do treat the numerical modeling of plasmonic structures as such, mainly the FDTD or the FEM technique is applied [37], [38]. Very few papers consider the IE-MoM technique, and if they do, it is the Surface IE technique [39], [40], [41]. In [42] Chremmos uses a magnetic type scalar integral equation to

to customize simulation. Data can be exported to Matlab or in ASCII.

manual meshing mode. A Graphical User Interface is available.

*4.2.2. Other non-commercial tools found in literature* 

occurring and occasionally just mention the numerical tool used.

**4.2. Non-commercial software tools** 

*4.2.1. MAGMAS 3D [31]: IE-MoM* 

found.

### *4.1.4. CST [30]: FIT*

CST Microwave Studio (CST MWS) is based on the finite integration technique (FIT). It allows to choose the time domain as well as the frequency domain approach. Despite the presence of transient, eigenmode, and frequency domain solvers within CST MWS, the transient solver is considered as the flag ship module. The Time Domain Solver calculates the broadband behavior of electromagnetic devices in one simulation run with an arbitrarily fine frequency resolution. The modeling of curved structures using the Perfect Boundary Approximation® technique and the modeling of thin perfectly electric conducting sheets with the Thin Sheet Technique® tries to cope with the typical difficulties inherent to FDTD methods for classical structures. Several mesh types can be chosen. The automatic mesh generator detects the important points inside the structure (fixpoints) and locates mesh nodes there. The user can manually add fixpoints on a structure, as well as fully control the number of mesh lines in each coordinate with regards to the specified wavelength. Energy based adaptation of the mesh allows to refine it in a predefined number of passes, providing a mesh refinement of sophisticated design features for the price of longer overall simulation times. Although the FIT in principle can handle material parameters changing over the dielectric volumes defined, this is not implemented yet.

CST, as a general purpose software package being a real competitor for HFSS in the traditional microwave field, has gained a lot of popularity in the last few years. Also for the analysis and design of plasmonic structures, more and more results obtained with CST can be found in literature (just google the words CST and plasmonics). A problem sometimes observed with CST is a ripple in the frequency response in case the tool settings are not appropriate. This is due to the fact that the flagship of CST is inherently a time domain solver.

### *4.1.5. Lumerical [15]: FDTD*

Lumerical is the leading software tool in the plasmonics and photonics community. It is based on the FDTD algorithm, and can be used for 2D and 3D topologies. The topology has to be generated in the GUI or can be imported with GDSII/SEM files. Basic shapes include triangles, rectangular blocks, cylinders, conic surfaces, polygons, rings, user-defined (parametric) surfaces, spheres and pyramids. Several boundary conditions can be used: absorbing (PML), periodic, Bloch, symmetric, asymmetric, and metal boundaries. A nonuniform mesh can be used and automesh algorithms are provided. Since it is a dedicated time domain solver, it has a sophisticated library of Lorentz, Drude, Debye models for the material parameters. Excitation of the structure is possible with waveguide sources, dipoles, plane waves, focused beams and diffraction-limited spots. A scripting language is available to customize simulation. Data can be exported to Matlab or in ASCII.

### **4.2. Non-commercial software tools**

### *4.2.1. MAGMAS 3D [31]: IE-MoM*

34 Plasmonics – Principles and Applications

*4.1.4. CST [30]: FIT* 

solver.

*4.1.5. Lumerical [15]: FDTD* 

decomposition techniques and fast solvers.

dielectric volumes defined, this is not implemented yet.

tools, propagation mode tools, and resonance mode tools. The time-harmonic problems can be formulated in 1D, 2D and 3D. The topologies can be isolated or may occur in periodic patterns, or a mixture of both. Interesting is that dedicated tools for problems posed on cylindrically symmetric geometries are available. The electromagnetic fields are discretized with higher order edge elements. JCMsuite contains an automatic mesh generator (with adaptive features), goal-oriented error estimators for adaptive grid refinement, domain-

CST Microwave Studio (CST MWS) is based on the finite integration technique (FIT). It allows to choose the time domain as well as the frequency domain approach. Despite the presence of transient, eigenmode, and frequency domain solvers within CST MWS, the transient solver is considered as the flag ship module. The Time Domain Solver calculates the broadband behavior of electromagnetic devices in one simulation run with an arbitrarily fine frequency resolution. The modeling of curved structures using the Perfect Boundary Approximation® technique and the modeling of thin perfectly electric conducting sheets with the Thin Sheet Technique® tries to cope with the typical difficulties inherent to FDTD methods for classical structures. Several mesh types can be chosen. The automatic mesh generator detects the important points inside the structure (fixpoints) and locates mesh nodes there. The user can manually add fixpoints on a structure, as well as fully control the number of mesh lines in each coordinate with regards to the specified wavelength. Energy based adaptation of the mesh allows to refine it in a predefined number of passes, providing a mesh refinement of sophisticated design features for the price of longer overall simulation times. Although the FIT in principle can handle material parameters changing over the

CST, as a general purpose software package being a real competitor for HFSS in the traditional microwave field, has gained a lot of popularity in the last few years. Also for the analysis and design of plasmonic structures, more and more results obtained with CST can be found in literature (just google the words CST and plasmonics). A problem sometimes observed with CST is a ripple in the frequency response in case the tool settings are not appropriate. This is due to the fact that the flagship of CST is inherently a time domain

Lumerical is the leading software tool in the plasmonics and photonics community. It is based on the FDTD algorithm, and can be used for 2D and 3D topologies. The topology has to be generated in the GUI or can be imported with GDSII/SEM files. Basic shapes include triangles, rectangular blocks, cylinders, conic surfaces, polygons, rings, user-defined (parametric) surfaces, spheres and pyramids. Several boundary conditions can be used: absorbing (PML), periodic, Bloch, symmetric, asymmetric, and metal boundaries. A nonuniform mesh can be used and automesh algorithms are provided. Since it is a dedicated MAGMAS 3D is the IE-MoM code developed at the Katholieke Universiteit Leuven, Belgium. It was originally developed in cooperation with the European Space Agency for planar and quasi-planar antenna and scattering structures embedded in a multilayered dielectric background medium and operating in the microwave frequency range. Starting from the topology considered, frequency band needed, and type of excitation, it calculates the network, radiation, and scattering characteristics of the structure under consideration. It is based on a full-wave Mixed-Potential formulation of Electric Field Integral Equations (MP-EFIE), originally applied only to 2D surface currents [22], [32]. New theoretical techniques were developed and implemented within the framework: the Expansion Wave Concept (EWC) [33], [34], the Dipole Modeling Technique (DMT) [35], special deembedding procedures [36], etc.. In 2007, it was extended with the capability to handle volumetric 3D currents with the VIE technique introduced in [21]. This was crucial in view of its application in plasmonics. In 2009, this led to the first verified simulated results in this field [25]. Now, it is extensively used for plasmonics. The MAGMAS mesh is based on a combination of rectangular and triangular mesh cells. Full mesh control is available in manual meshing mode. A Graphical User Interface is available.

To the best knowledge of the author, to date (June 2012), it is the only IE-MoM framework able to handle arbitrary plasmonic structures embedded in a multilayered environment, with the very general and flexible VIE technique. No similar commercial solvers exist, in contrast to the situation at microwave frequencies, where a multitude of IE solvers can be found.

#### *4.2.2. Other non-commercial tools found in literature*

At this moment plasmonics is a mainly experimentally driven research field. Whereas for example in the field of traditional antenna research, numerous papers can be found on modeling as such, this is much less the case in plasmonics, especially in the high impact journals. Most papers are concerned with the description of the physical phenomena occurring and occasionally just mention the numerical tool used.

Most researchers in the plasmonics field use the commercial tools available, and as far as I can see, the main one is Lumerical. In the publications that do treat the numerical modeling of plasmonic structures as such, mainly the FDTD or the FEM technique is applied [37], [38]. Very few papers consider the IE-MoM technique, and if they do, it is the Surface IE technique [39], [40], [41]. In [42] Chremmos uses a magnetic type scalar integral equation to

#### 36 Plasmonics – Principles and Applications

describe surface plasmon scattering by rectangular dielectric channel discontinuities. Even the use of the magnetic current formalism to describe holes, classical at microwave frequencies, already has been used in plasmonics [43]. However, these dedicated developments cannot be categorized under the title "analysis framework" in the sense that they do not allow to handle a wide range of different topologies.

Computational Electromagnetics in Plasmonics 37

and FDTD can reach a high level of accuracy only with a high discretization. In 2D, this is readily affordable on present-day computer systems. There, these techniques have a clear advantage in terms of speed, matrix size, and accuracy. In 3D, due to the rocketing size of the problem, this may not always be that obvious. The advantages may thus disappear in a full 3D analysis of geometrically complicated structures. This happens because matrices become denser or more ill-conditioned. The main conclusion of the paper is that the most efficient method depends on the problem dimension and complexity. This is a similar conclusion as was reached in [46] within the context of the analysis and design of planar

**5.2. Comparison between differential and integral equation techniques** 

In this section, the differential and integral equation techniques are compared. This is done by choosing a representative solver from each category and using it in the analysis of a basic plasmonic topology. In the category differential techniques, Lumerical is chosen, as it is the most widespread commercial solver in the plasmonic community. Since, as far as we know, there are no commercial integral equation solvers in the plasmonics area, MAGMAS 3D is chosen in this category, the in-house developed solver at Katholieke Universiteit Leuven. The basic topology selected is the plasmonic dipole. It is a structure that can function as a scatterer, but also as a real nano-antenna. The main result of this section is that it draws conclusions on different computational aspects involved in the modeling of nanostructures

The plasmonic dipole is depicted in Fig. 3. It is a structure consisting of two nanorods with a gap in between. When the gap is "shortcircuited", it functions as a single rod, when the gap is open, it functions as a device generating an enhanced electric field there. When it is connected to other nanocircuits, it may function as a nano-antenna, both in transmit (Fig. 3b), and in receive (Fig. 3c). The two rods may be fabricated from metals like Au, Ag, Cu, Al, Cr, etc.. In this section the width W and height H of the dipole are kept constant at 40 nm. This value is well-chosen, since it is a value which can be fabricated with sufficient accuracy using present-day nanofabrication technology. The gap width can have different values, depending on the case considered. Note however that 10 nm is about the minimum that can

It cannot be over-emphasized that Computational Electromagnetics solvers in general produce correct results consistently for a multitude of different structures only if two

1. The user has sufficient general background in the field, and a more specific knowledge about the solution technique implemented in the solver that he is using. If not, there is a huge danger that the solver is not used in a proper way. For example, although this is not necessary to "operate" the solver, if the user has no basic knowledge and does not

antennas.

with the two solvers.

*5.2.1. Plasmonic dipole* 

*5.2.1.1. Quality of the input* 

conditions are met:

be fabricated with reasonable accuracy nowadays.

### **5. Benchmarking**

Nowadays, physicists and engineers rely heavily on highly specialized full wave electromagnetic field solvers to analyze, develop, and optimize their designs. Computeraided analysis and optimization have replaced the design process of iterative experimental modifications of the initial design. It is evident that the underlying solution method for a software tool may significantly influence the efficiency and accuracy by which certain structure types are analyzed. Nevertheless, the commercial focus increasingly switches from such key theoretical considerations to improvements in the area of layout tools and systemlevel design tools. Therefore, users may get the wrong impression that a given solver is automatically suited to solve any kind of problem with arbitrary precision. This is of course not true.

This section verifies the plausibility of such expectations by presenting a benchmark study for a few plasmonic structures. The study focuses on the capabilities and limitations of the applied EM modeling techniques that usually remain hidden for the user.

### **5.1. Benchmarking in literature**

In literature, not many comparisons between solvers in the field of plasmonics can be found. In [44] Hoffmann et al. consider a single plasmonic topology, a pair of Au spheres, and analyze this structure with several electromagnetic field solvers: COMSOL Multiphysics (FEM), JCMsuite (FEM), HFSS (FEM), and CST Microwave Studio (the FEM solver available in the Studio is used, not the flagship FDTD solver). The output parameter considered is the electric field strength in between the two spheres. Note that all solvers tested are based on the Finite Elements technique in the frequency domain, so it can be expected that the main issues with these solvers for this topology are the same. Very interesting is the fact that in the paper itself CST is categorized as "inaccurate". However, afterwards this was corrected. It seems that wrong material settings were used, and after correction accurate results were obtained. This issue led to the fact that this benchmark example is now presented as a reference example on the CST website.

In [45], several numerical methods are tested for 2D plasmonic nanowire structures: not only the Finite Element Method (FEM) and the Finite Difference Time-Domain (FDTD) technique, but also less "commercial" methods like the Multiple Multipole Program (MMP), the Method of Auxiliary Sources (MAS), and the Mesh-less Boundary Integral Equation (BIE) method are tested. By comparing the results, several conclusions can be drawn about their applicability and accuracy for plasmonic topologies. Differential techniques like FEM and FDTD can reach a high level of accuracy only with a high discretization. In 2D, this is readily affordable on present-day computer systems. There, these techniques have a clear advantage in terms of speed, matrix size, and accuracy. In 3D, due to the rocketing size of the problem, this may not always be that obvious. The advantages may thus disappear in a full 3D analysis of geometrically complicated structures. This happens because matrices become denser or more ill-conditioned. The main conclusion of the paper is that the most efficient method depends on the problem dimension and complexity. This is a similar conclusion as was reached in [46] within the context of the analysis and design of planar antennas.

### **5.2. Comparison between differential and integral equation techniques**

In this section, the differential and integral equation techniques are compared. This is done by choosing a representative solver from each category and using it in the analysis of a basic plasmonic topology. In the category differential techniques, Lumerical is chosen, as it is the most widespread commercial solver in the plasmonic community. Since, as far as we know, there are no commercial integral equation solvers in the plasmonics area, MAGMAS 3D is chosen in this category, the in-house developed solver at Katholieke Universiteit Leuven. The basic topology selected is the plasmonic dipole. It is a structure that can function as a scatterer, but also as a real nano-antenna. The main result of this section is that it draws conclusions on different computational aspects involved in the modeling of nanostructures with the two solvers.

### *5.2.1. Plasmonic dipole*

36 Plasmonics – Principles and Applications

**5. Benchmarking** 

not true.

**5.1. Benchmarking in literature** 

reference example on the CST website.

describe surface plasmon scattering by rectangular dielectric channel discontinuities. Even the use of the magnetic current formalism to describe holes, classical at microwave frequencies, already has been used in plasmonics [43]. However, these dedicated developments cannot be categorized under the title "analysis framework" in the sense that

Nowadays, physicists and engineers rely heavily on highly specialized full wave electromagnetic field solvers to analyze, develop, and optimize their designs. Computeraided analysis and optimization have replaced the design process of iterative experimental modifications of the initial design. It is evident that the underlying solution method for a software tool may significantly influence the efficiency and accuracy by which certain structure types are analyzed. Nevertheless, the commercial focus increasingly switches from such key theoretical considerations to improvements in the area of layout tools and systemlevel design tools. Therefore, users may get the wrong impression that a given solver is automatically suited to solve any kind of problem with arbitrary precision. This is of course

This section verifies the plausibility of such expectations by presenting a benchmark study for a few plasmonic structures. The study focuses on the capabilities and limitations of the

In literature, not many comparisons between solvers in the field of plasmonics can be found. In [44] Hoffmann et al. consider a single plasmonic topology, a pair of Au spheres, and analyze this structure with several electromagnetic field solvers: COMSOL Multiphysics (FEM), JCMsuite (FEM), HFSS (FEM), and CST Microwave Studio (the FEM solver available in the Studio is used, not the flagship FDTD solver). The output parameter considered is the electric field strength in between the two spheres. Note that all solvers tested are based on the Finite Elements technique in the frequency domain, so it can be expected that the main issues with these solvers for this topology are the same. Very interesting is the fact that in the paper itself CST is categorized as "inaccurate". However, afterwards this was corrected. It seems that wrong material settings were used, and after correction accurate results were obtained. This issue led to the fact that this benchmark example is now presented as a

In [45], several numerical methods are tested for 2D plasmonic nanowire structures: not only the Finite Element Method (FEM) and the Finite Difference Time-Domain (FDTD) technique, but also less "commercial" methods like the Multiple Multipole Program (MMP), the Method of Auxiliary Sources (MAS), and the Mesh-less Boundary Integral Equation (BIE) method are tested. By comparing the results, several conclusions can be drawn about their applicability and accuracy for plasmonic topologies. Differential techniques like FEM

applied EM modeling techniques that usually remain hidden for the user.

they do not allow to handle a wide range of different topologies.

The plasmonic dipole is depicted in Fig. 3. It is a structure consisting of two nanorods with a gap in between. When the gap is "shortcircuited", it functions as a single rod, when the gap is open, it functions as a device generating an enhanced electric field there. When it is connected to other nanocircuits, it may function as a nano-antenna, both in transmit (Fig. 3b), and in receive (Fig. 3c). The two rods may be fabricated from metals like Au, Ag, Cu, Al, Cr, etc.. In this section the width W and height H of the dipole are kept constant at 40 nm. This value is well-chosen, since it is a value which can be fabricated with sufficient accuracy using present-day nanofabrication technology. The gap width can have different values, depending on the case considered. Note however that 10 nm is about the minimum that can be fabricated with reasonable accuracy nowadays.

#### *5.2.1.1. Quality of the input*

It cannot be over-emphasized that Computational Electromagnetics solvers in general produce correct results consistently for a multitude of different structures only if two conditions are met:

1. The user has sufficient general background in the field, and a more specific knowledge about the solution technique implemented in the solver that he is using. If not, there is a huge danger that the solver is not used in a proper way. For example, although this is not necessary to "operate" the solver, if the user has no basic knowledge and does not realize how the discretization / meshing scheme works, it is impossible to understand and assess the effect of a proper meshing. In the majority of the cases, there will be no problem with that, and automatic meshing schemes will be able to produce good results. However, in unusual and/or challenging cases, often met in scientific research, this issue will be a crucial factor.

Computational Electromagnetics in Plasmonics 39

seen that there is an excellent agreement in the trend of the curves, i.e. the change of input impedance with frequency for different widths of the source current filament in MAGMAS (5 nm and 13.33 nm) is the same as the one in Lumerical. However, in general there is an offset. Only for a specific width of the current filament in MAGMAS, nl. 5 nm, the two excitation mechanisms become almost completely equivalent and even the offset between the two solvers disappears. A similar agreement was also observed for both gold and silver

The main conclusion of this section is that there is still an issue with the modeling of (localized) feeds, especially in Lumerical. The reason is that the plasmonics community is used to consider objects as scatterers, excited by an incident wave. In this case, the models implemented are satisfactory. However, in view of designing nanocircuits, localized feeds and impedances derived from them become important. The basic rule is that the feeding model has to correspond as closely as possible to the actual feed topology that is going to be used later in practice. In this respect, the localized feeding models available in Lumerical

**Figure 4.** Comparison of impedances obtained with MAGMAS and Lumerical for an aluminum dipole on a glass substrate, L = 200 nm and gap = 10 nm. Two widths *w* are used in MAGMAS for the source

<sup>300</sup> <sup>400</sup> <sup>500</sup> <sup>600</sup> <sup>700</sup> <sup>800</sup> <sup>900</sup> <sup>1000</sup> <sup>1100</sup> <sup>1200</sup> �150

Wavelength, nm

Lumerical

MAGMAS:W=5nm MAGMAS:W=13.33nm

<sup>300</sup> <sup>400</sup> <sup>500</sup> <sup>600</sup> <sup>700</sup> <sup>800</sup> <sup>900</sup> <sup>1000</sup> <sup>1100</sup> <sup>1200</sup> <sup>0</sup>

MAGMAS:W=5nm MAGMAS:W=13.33nm Lumerical

Wavelength, nm

current filament. R and X stand for real and imaginary part.

�100

�50

0

X, Ω

50

100

150

R, Ω

dipoles.

today are unsatisfactory.

2. The user needs to have sufficient knowledge about the actual tool that he is using and its peculiarities, more specific about the way a certain problem has to be imported into the solver. A very good illustration of this has already been addressed in section 5.1 where CST was actually not properly used and produced incorrect results that were published in open literature. Afterwards, this was corrected and the problem disappeared. However, following the proper line of reasoning is not always straightforward. Even experienced researchers sometimes easily can make mistakes. It is the opinion of the author that this is happening too much, even in a lot of peer-reviewed scientific papers.

**Figure 3.** The dipole model studied. a) W and H are set equal to 40 nm and the gap G can vary. b) The dipole as transmitting antenna with a model for the feeding structure located in the gap. c) The dipole as receiving antenna excited by a plane wave.

#### *5.2.1.2. Modeling of the feed*

The first parameter studied is the impedance that is seen when the nanodevice operates as a transmitting or receiving antenna. In Fig. 4, the impedance simulated with both MAGMAS and Lumerical is given for an Al dipole of 200 nm length and with a gap of 10 nm. The input impedance is calculated as *Z = V/I*, where *Z*, *V* and *I* are the input impedance, the voltage over the gap G, and the source current, respectively. The exciting source is located in the middle of the gap. MAGMAS uses a physical current filament of finite width *w*. Lumerical does not have a built-in physically feasible current source model. It uses a (less realistic) pure electrical dipole source, which means that the current *I* has to be evaluated "manually" by the user based on the magnetic field distribution and Ampere's law. This also means that the impedance will be somewhat depending on the resolution chosen during meshing. It is seen that there is an excellent agreement in the trend of the curves, i.e. the change of input impedance with frequency for different widths of the source current filament in MAGMAS (5 nm and 13.33 nm) is the same as the one in Lumerical. However, in general there is an offset. Only for a specific width of the current filament in MAGMAS, nl. 5 nm, the two excitation mechanisms become almost completely equivalent and even the offset between the two solvers disappears. A similar agreement was also observed for both gold and silver dipoles.

38 Plasmonics – Principles and Applications

this issue will be a crucial factor.

as receiving antenna excited by a plane wave.

*5.2.1.2. Modeling of the feed* 

realize how the discretization / meshing scheme works, it is impossible to understand and assess the effect of a proper meshing. In the majority of the cases, there will be no problem with that, and automatic meshing schemes will be able to produce good results. However, in unusual and/or challenging cases, often met in scientific research,

2. The user needs to have sufficient knowledge about the actual tool that he is using and its peculiarities, more specific about the way a certain problem has to be imported into the solver. A very good illustration of this has already been addressed in section 5.1 where CST was actually not properly used and produced incorrect results that were published in open literature. Afterwards, this was corrected and the problem disappeared. However, following the proper line of reasoning is not always straightforward. Even experienced researchers sometimes easily can make mistakes. It is the opinion of the author that this is

**Figure 3.** The dipole model studied. a) W and H are set equal to 40 nm and the gap G can vary. b) The dipole as transmitting antenna with a model for the feeding structure located in the gap. c) The dipole

The first parameter studied is the impedance that is seen when the nanodevice operates as a transmitting or receiving antenna. In Fig. 4, the impedance simulated with both MAGMAS and Lumerical is given for an Al dipole of 200 nm length and with a gap of 10 nm. The input impedance is calculated as *Z = V/I*, where *Z*, *V* and *I* are the input impedance, the voltage over the gap G, and the source current, respectively. The exciting source is located in the middle of the gap. MAGMAS uses a physical current filament of finite width *w*. Lumerical does not have a built-in physically feasible current source model. It uses a (less realistic) pure electrical dipole source, which means that the current *I* has to be evaluated "manually" by the user based on the magnetic field distribution and Ampere's law. This also means that the impedance will be somewhat depending on the resolution chosen during meshing. It is

happening too much, even in a lot of peer-reviewed scientific papers.

The main conclusion of this section is that there is still an issue with the modeling of (localized) feeds, especially in Lumerical. The reason is that the plasmonics community is used to consider objects as scatterers, excited by an incident wave. In this case, the models implemented are satisfactory. However, in view of designing nanocircuits, localized feeds and impedances derived from them become important. The basic rule is that the feeding model has to correspond as closely as possible to the actual feed topology that is going to be used later in practice. In this respect, the localized feeding models available in Lumerical today are unsatisfactory.

**Figure 4.** Comparison of impedances obtained with MAGMAS and Lumerical for an aluminum dipole on a glass substrate, L = 200 nm and gap = 10 nm. Two widths *w* are used in MAGMAS for the source current filament. R and X stand for real and imaginary part.

#### 40 Plasmonics – Principles and Applications

#### *5.2.1.3. Modeling of losses*

The importance of taking into account losses is obvious. Metals at plasmonic frequencies may be very lossy and neglecting this would result in totally erroneous results. However, there is a difference between time domain and frequency domain solvers. Since Lumerical is an FDTD based tool, it works in time domain. This means that it cannot work directly with (complex) permittivities and permeabilities. The material data have to be transformed from functions of frequency into functions of time. This is a very difficult process, involving a convolution, and is subject to various constraints. Local curve fitting of permittivities is very difficult to implement in a time domain solver. Lumerical fits the sampled measurement data with its own multi-coefficient material model, which "provides a superior fit compared to the standard material models (such as Lorentz, Debye, Drude, etc )". However, this produces artefacts, as is illustrated in Fig. 5. There, the relative error is given between the measured imaginary part of the permittivity (which represents losses) and the value obtained for this imaginary part after fitting, for three metals. It is seen that for gold and aluminum, the error is in the order of 20 %, which is already quite elevated. However, in the case of silver, the error reaches a full 100 %, which is unacceptable. In some cases, it will result in a totally wrong prediction for example of the radiation efficiency of nano-antennas [47]. If one forces the material characteristics to be exactly the same in the two solvers, for example by using the interpolated values also in MAGMAS, instead of the real measured ones, the two solvers produce very similar results, see Fig. 6.

Computational Electromagnetics in Plasmonics 41

Lumerical FDTD Solver MAGMAS MOM Solver

**Figure 6.** Comparison between FDTD and IE-MoM for a gold dipole in free space, L = 250 nm and the

400 600 800 1000 1200 1400

Wavelength,nm

It is well known that the mesh quality and resolution are key factors in the accuracy of any solver. The electromagnetic coupling between nearby segments may differ considerably due to the specific meshing used, especially in parts of the structure where rapid topological changes occur. The question is whether the meshes used in the solvers are adequate. This issue was investigated by performing a convergence study in terms of the resolution of the

The extinction cross section of a gold monomer (no gap, or in other words G = 0 nm) calculated using MAGMAS 3D with different meshes is plotted in Fig. 7(a). The analysis of these data demonstrates clearly the stability of IE-MoM. The results obtained even with a very rough mesh 6x1x1 (41.7 *nm* x 40 *nm* x 40 *nm*) provide already a very good estimation of the antenna resonance properties. The differences in extinction cross section calculated with the rough and fine meshes are almost negligibly small. This result is obtained thanks to the variational stability of IE-MoM. On the adjacent figure Fig. 7(b) extinction cross sections are calculated using Lumerical. As expected, in general the results obtained with the FDTD method depend considerably on the chosen mesh. The calculated wavelengths as a function of the mesh cell size are plotted in Fig. 7(c). It should be also noted that in contrast to Fig. 7(a) (IE-MoM) in Fig. 7(b) (FDTD) not only peak positions but also their levels depend

The conclusion is that in Lumerical a fine mesh is mandatory for reliable calculations of the monomer resonant wavelength. Thanks to its variational approach, in IE-MoM, a highly

dense mesh is not always needed, especially when scattering problems are studied.

gap is 10 nm. Identical material parameters are used in both solvers.

0

0.2

0.4

0.6

Radiation Efficiency

0.8

1

*5.2.1.4. Meshing* 

meshes [48].

clearly on the mesh.

The main conclusion is that a highly accurate prediction of losses in Lumerical is still an issue. The material models have to be further refined, especially for silver.

**Figure 5.** Relative error between interpolated and measured imaginary part of permittivity for silver, gold and aluminum.

**Figure 6.** Comparison between FDTD and IE-MoM for a gold dipole in free space, L = 250 nm and the gap is 10 nm. Identical material parameters are used in both solvers.

#### *5.2.1.4. Meshing*

40 Plasmonics – Principles and Applications

The importance of taking into account losses is obvious. Metals at plasmonic frequencies may be very lossy and neglecting this would result in totally erroneous results. However, there is a difference between time domain and frequency domain solvers. Since Lumerical is an FDTD based tool, it works in time domain. This means that it cannot work directly with (complex) permittivities and permeabilities. The material data have to be transformed from functions of frequency into functions of time. This is a very difficult process, involving a convolution, and is subject to various constraints. Local curve fitting of permittivities is very difficult to implement in a time domain solver. Lumerical fits the sampled measurement data with its own multi-coefficient material model, which "provides a superior fit compared to the standard material models (such as Lorentz, Debye, Drude, etc )". However, this produces artefacts, as is illustrated in Fig. 5. There, the relative error is given between the measured imaginary part of the permittivity (which represents losses) and the value obtained for this imaginary part after fitting, for three metals. It is seen that for gold and aluminum, the error is in the order of 20 %, which is already quite elevated. However, in the case of silver, the error reaches a full 100 %, which is unacceptable. In some cases, it will result in a totally wrong prediction for example of the radiation efficiency of nano-antennas [47]. If one forces the material characteristics to be exactly the same in the two solvers, for example by using the interpolated values also in MAGMAS, instead of the real measured

The main conclusion is that a highly accurate prediction of losses in Lumerical is still an

Silver Gold

Aluminum

**Figure 5.** Relative error between interpolated and measured imaginary part of permittivity for silver,

200 300 400 500 600 700 800 900 1000

Frequency, THz

ones, the two solvers produce very similar results, see Fig. 6.

issue. The material models have to be further refined, especially for silver.

*5.2.1.3. Modeling of losses* 

gold and aluminum.

�1 �0.8 �0.6 �0.4 �0.2 0 0.2 0.4 0.6 0.8 1

Relative Error

It is well known that the mesh quality and resolution are key factors in the accuracy of any solver. The electromagnetic coupling between nearby segments may differ considerably due to the specific meshing used, especially in parts of the structure where rapid topological changes occur. The question is whether the meshes used in the solvers are adequate. This issue was investigated by performing a convergence study in terms of the resolution of the meshes [48].

The extinction cross section of a gold monomer (no gap, or in other words G = 0 nm) calculated using MAGMAS 3D with different meshes is plotted in Fig. 7(a). The analysis of these data demonstrates clearly the stability of IE-MoM. The results obtained even with a very rough mesh 6x1x1 (41.7 *nm* x 40 *nm* x 40 *nm*) provide already a very good estimation of the antenna resonance properties. The differences in extinction cross section calculated with the rough and fine meshes are almost negligibly small. This result is obtained thanks to the variational stability of IE-MoM. On the adjacent figure Fig. 7(b) extinction cross sections are calculated using Lumerical. As expected, in general the results obtained with the FDTD method depend considerably on the chosen mesh. The calculated wavelengths as a function of the mesh cell size are plotted in Fig. 7(c). It should be also noted that in contrast to Fig. 7(a) (IE-MoM) in Fig. 7(b) (FDTD) not only peak positions but also their levels depend clearly on the mesh.

The conclusion is that in Lumerical a fine mesh is mandatory for reliable calculations of the monomer resonant wavelength. Thanks to its variational approach, in IE-MoM, a highly dense mesh is not always needed, especially when scattering problems are studied.

Computational Electromagnetics in Plasmonics 43

(10 processors in parallel)

IE-MoM FDTD

Similarly to a pebble hitting a water surface, the energy deposition by femtosecond laser pulses on a gold surface can produce local melting and back-jet. Contrary to water though, a gold surface can be nanopatterned and the excitation of surface plasmon resonances leads to the appearance of hotspots, literally. This was explicitly proven for the first time in a nanopatterned gold surface, composed of a G shaped periodic structure. It was observed there that laser-induced melting and back-jet occur precisely in the plasmonic hotspots. The aid of computational electromagnetics to rigorously analyze the structure by predicting the hotspots was crucial in explaining this phenomenon. Much more details can be found in [49]. It is clear that this type of insight into the basic interaction mechanisms of light with

Another example of what can be reached by applying computational electromagnetics in

**Figure 8.** Illuminating the sample with femtosecond laser pulses produces a polarization dependent pattern of nanobumps with sharp tips on the sample surface [49]. Left: horizontally polarized light,

6 x 6 x 1 5.4 6 x 1 x 1 38 12 x 2 x 2 10.4 12 x 2 x 2 55 25 x 4 x 4 33.6 25 x 4 x 4 225 50 x 8 x 6 238 50 x 8 x 8 672

frequency point (s) Mesh Total calculation time (s)

100 x 16 x 16 5893 250 x 40 x 40 20100

Mesh Calculation time per

MoM: Intel(R) Xeon(R) CPU E5335 @ 2 GHz, 32 Gb memory FDTD: d1585g6 4x hc CPU @ 2.8 GHz, 128 Gb memory

plasmonics can be found in [50].

right: vertically polarized light.

**Table 1.** Calculation times as a function of the mesh for the gold monomer.

matter at the nanoscale opens up new possibilities for applications.

**Figure 7.** a and b: extinction cross section of a 250 nm x 40 nm x 40 nm gold monomer on a substrate with n=1.5 (a. MAGMAS 3D, b. Lumerical), c. convergence performance of both solvers.

#### *5.2.1.5. Calculation speed*

Information on the calculation times for the convergence study of the previous section is given in Table 1. It has to be emphasized that calculation times cannot be directly compared. The computers used were different, and the FDTD tool used 10 processors in parallel. Further, since it is a time domain technique, a calculation time per frequency point cannot be given for FDTD. Nevertheless, combining the data in this table with the data concerning the convergence, given in Fig. 7, it is easily seen that IE-MoM outperforms FDTD for this particular case.

#### *5.2.2. Plasmonic nanojets*

In this section it is illustrated what can be reached with computational electromagnetics in the field of plasmonics. The results in this section are extracted from the paper [49], where full details can be found.


MoM: Intel(R) Xeon(R) CPU E5335 @ 2 GHz, 32 Gb memory FDTD: d1585g6 4x hc CPU @ 2.8 GHz, 128 Gb memory

42 Plasmonics – Principles and Applications

(a) (b)

1200

1250

1300

Resonant wavelength [nm]

1350

*5.2.1.5. Calculation speed* 

particular case.

*5.2.2. Plasmonic nanojets* 

full details can be found.

**Figure 7.** a and b: extinction cross section of a 250 nm x 40 nm x 40 nm gold monomer on a substrate

0 10 20 30 40

Mesh size along polarization [nm]

Information on the calculation times for the convergence study of the previous section is given in Table 1. It has to be emphasized that calculation times cannot be directly compared. The computers used were different, and the FDTD tool used 10 processors in parallel. Further, since it is a time domain technique, a calculation time per frequency point cannot be given for FDTD. Nevertheless, combining the data in this table with the data concerning the convergence, given in Fig. 7, it is easily seen that IE-MoM outperforms FDTD for this

In this section it is illustrated what can be reached with computational electromagnetics in the field of plasmonics. The results in this section are extracted from the paper [49], where

with n=1.5 (a. MAGMAS 3D, b. Lumerical), c. convergence performance of both solvers.

MAGMAS3D Lumerical

**Table 1.** Calculation times as a function of the mesh for the gold monomer.

Similarly to a pebble hitting a water surface, the energy deposition by femtosecond laser pulses on a gold surface can produce local melting and back-jet. Contrary to water though, a gold surface can be nanopatterned and the excitation of surface plasmon resonances leads to the appearance of hotspots, literally. This was explicitly proven for the first time in a nanopatterned gold surface, composed of a G shaped periodic structure. It was observed there that laser-induced melting and back-jet occur precisely in the plasmonic hotspots. The aid of computational electromagnetics to rigorously analyze the structure by predicting the hotspots was crucial in explaining this phenomenon. Much more details can be found in [49]. It is clear that this type of insight into the basic interaction mechanisms of light with matter at the nanoscale opens up new possibilities for applications.

Another example of what can be reached by applying computational electromagnetics in plasmonics can be found in [50].

**Figure 8.** Illuminating the sample with femtosecond laser pulses produces a polarization dependent pattern of nanobumps with sharp tips on the sample surface [49]. Left: horizontally polarized light, right: vertically polarized light.

Computational Electromagnetics in Plasmonics 45

I would like to end this chapter with the following guideline. The use of two different solvers, based on different theoretical methods (integral and differential) may provide an excellent means to characterize the quality of simulation results. If the two results are in good agreement, it is highly likely that the results are correct. If the two results are in disagreement, a deeper investigation of the structure and its modeling is absolutely

The author gratefully acknowledges the following persons: Zhongkun Ma for the impedance and radiation efficiency comparison between MAGMAS and Lumerical, Dr. Vladimir Volski for the convergence study with MAGMAS and Lumerical, and V. K. Valev

We also would like to express our gratitude to the KU Leuven Methusalem project and its beneficiary Prof. V. V. Moshchalkov, whose activities created a real stimulus to become personally active in this field, and the fund for scientific research Flanders (FWO-V) for the

[1] I. Ahmed, E. H. Khoo, E. Li, and R. Mittra, "A hybrid approach for solving coupled Maxwell and Schrödinger equations arising in the simulation of nano-devices", IEEE

[2] D. B. Davidson,"A review of important recent developments in full-wave CEM for RF and microwave engineering," IEEE 3rd Int. Conf. Comp. Electromagnetics and Its

[3] C. W. Townbridge and J. K. Sykulski,"Some Key Developments in Computational Electromagnetics and Their Attribution," IEEE Antennas and Propag. Magazine, vol.

[4] E. K. Miller,"A Selective Survey of Computational Electromagnetics," IEEE Trans

[5] F. Peterson, S. L. Ray, and R. Mittra, "Computational methods for electromagnetics",

[6] Awadhiya, P. Barba, and L. Kempel, "Finite-element method programming made easy???", IEEE Antennas and Propagation Magazine, vol. 45, pp. 73 – 79, Aug. 2003. [7] J. M. Jin, The Finite Element Method in Electromagnetics, second edition, John Willey &

Antennas and Wireless Propagation Letters, Vol. 9, pp. 914- 917, 2010.

Antennas Propag, vol. 36, no. 9, Sept. 1988, pp. 1281 – 1305

necessary.

**Author details** 

Guy A. E. Vandenbosch

**Acknowledgement** 

financial support.

**7. References** 

*Katholieke Universiteit Leuven, Belgium* 

for making available the plasmonic nanojet figures.

Applications, pp. PS/1-PS/4, Nov. 2004.

IEEE Press – Oxford University Press, 1998.

42, no. 6, pp. 503 – 508, Apr. 2006

Sons, Inc., New York, 2002.

**Figure 9.** The pattern of melted nanobumps matches that of the plasmonic currents, where Ohmic losses locally increase the temperature [49]. The electric field of light incident on the nanostructures causes the surface charges to oscillate in response to the direction of light polarization: horizontal, in (a), and vertical, in (b). The locations of the resulting electric currents are shown in (c) and (d) respectively. The current maxima are indicated with colored circles matching the pattern of nanobumps in Fig. 8.

### **6. Conclusions**

In this chapter, a brief introductory overview has been given on the use of existing computational techniques and software tools for the analysis and design of plasmonic structures. This field is booming, and many modeling techniques developed at lower frequencies, i.e. in the microwave range, are now being transferred to the plasmonics community, of course with the necessary adaptations. Although the differential equation techniques are the most widespread, both in scientific literature, and in the commercial scene, it is proven that integral equation techniques are a valid alternative. In some cases they are even superior. The author sees a lot of opportunities in this area.

In most cases, the tools are used to analyze structures. Few papers really make dedicated designs. However, this will change in the future. Inevitably, the plasmonics community will follow a similar patch as the traditional antenna community. Whereas 30 years ago, the design of an antenna was based on the accumulated expertise and know-how of many years, nowadays most antennas are designed almost in a single pass by experienced designers using commercial tools. This is possible since the tools in this field have sufficient accuracy and matureness in order to be able to do that. It is my belief that as soon as plasmonics will make the unavoidable shift from the physics to the engineering community, this last community will aim at conceiving applications which will involve designs heavily based on computational electromagnetics.

I would like to end this chapter with the following guideline. The use of two different solvers, based on different theoretical methods (integral and differential) may provide an excellent means to characterize the quality of simulation results. If the two results are in good agreement, it is highly likely that the results are correct. If the two results are in disagreement, a deeper investigation of the structure and its modeling is absolutely necessary.

### **Author details**

44 Plasmonics – Principles and Applications

**6. Conclusions** 

**Figure 9.** The pattern of melted nanobumps matches that of the plasmonic currents, where Ohmic losses locally increase the temperature [49]. The electric field of light incident on the nanostructures causes the surface charges to oscillate in response to the direction of light polarization: horizontal, in (a), and vertical, in (b). The locations of the resulting electric currents are shown in (c) and (d) respectively. The current maxima are indicated with colored circles matching the pattern of nanobumps in Fig. 8.

In this chapter, a brief introductory overview has been given on the use of existing computational techniques and software tools for the analysis and design of plasmonic structures. This field is booming, and many modeling techniques developed at lower frequencies, i.e. in the microwave range, are now being transferred to the plasmonics community, of course with the necessary adaptations. Although the differential equation techniques are the most widespread, both in scientific literature, and in the commercial scene, it is proven that integral equation techniques are a valid alternative. In some cases

In most cases, the tools are used to analyze structures. Few papers really make dedicated designs. However, this will change in the future. Inevitably, the plasmonics community will follow a similar patch as the traditional antenna community. Whereas 30 years ago, the design of an antenna was based on the accumulated expertise and know-how of many years, nowadays most antennas are designed almost in a single pass by experienced designers using commercial tools. This is possible since the tools in this field have sufficient accuracy and matureness in order to be able to do that. It is my belief that as soon as plasmonics will make the unavoidable shift from the physics to the engineering community, this last community will aim at conceiving applications which will involve designs heavily

they are even superior. The author sees a lot of opportunities in this area.

based on computational electromagnetics.

Guy A. E. Vandenbosch *Katholieke Universiteit Leuven, Belgium* 

### **Acknowledgement**

The author gratefully acknowledges the following persons: Zhongkun Ma for the impedance and radiation efficiency comparison between MAGMAS and Lumerical, Dr. Vladimir Volski for the convergence study with MAGMAS and Lumerical, and V. K. Valev for making available the plasmonic nanojet figures.

We also would like to express our gratitude to the KU Leuven Methusalem project and its beneficiary Prof. V. V. Moshchalkov, whose activities created a real stimulus to become personally active in this field, and the fund for scientific research Flanders (FWO-V) for the financial support.

### **7. References**


[8] J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite element method for electromagnetics, IEEE Press, Oxford University Press, 1997.

Computational Electromagnetics in Plasmonics 47

[25] G. A. E. Vandenbosch, V. Volski, N. Verellen, and V. V. Moshchalkov, "On the use of the method of moments in plasmonic applications", Radio Science, 46, RS0E02,

[26] R. Qiang, R. L. Chen, and J. Chen, "Modeling electrical properties of gold films at infrared frequency using FDTD method", *International Journal of Infrared and Millimeter* 

[30] CST GmbH – Computer Simulation Technology, CST Microwave Studio 2006 user

[31] G. A. E. Vandenbosch, MAGMAS 3D, http://www.esat.kuleuven.be/telemic/antennas/

[32] G. A. E. Vandenbosch and A. R. Van de Capelle, Mixed-potential integral expression formulation of the electric field in a stratified dielectric medium - application to the case of a probe current source, *IEEE Trans. Antennas Propagat*., vol. 40, pp. 806-817, July 1992. [33] F. J. Demuynck, G. A. E. Vandenbosch and A. R. Van de Capelle, The expansion wave concept, part I: efficient calculation of spatial Green's functions in a stratified dielectric

[34] G. A. E. Vandenbosch and F. J. Demuynck, The expansion wave concept, part II: a new way to model mutual coupling in microstrip antennas, *IEEE Trans. Antennas Propagat.*,

[35] B. L. A. Van Thielen and G. A. E. Vandenbosch, Method for the calculation of mutual coupling between discontinuities in Planar Circuits, *IEEE Trans. Microwave Theory and* 

[36] E. A. Soliman, G. A. E. Vandenbosch, E. Beyne and R.P. Mertens, Multimodal Characterization of Planar Microwave Structures, *IEEE Trans. Microwave Theory and* 

[37] R. Kappeler, D. Erni, X.-D. Cui, L, Novotny, "Field computations of optical antennas", *Journal of Computational and Theoretical Nanoscience*, vol. 4, no. 3, pp. 686-691, 2007. [38] R. Qiang, R. L. Chen, and J. Chen, "Modeling electrical properties of gold films at infrared frequency using FDTD method", *International Journal of Infrared and Millimeter* 

[39] B. Gallinet and O. J. Martin, "Scattering on plasmonic nanostructures arrays modeled with a surface integral formulation", *Photonics and Nanostructures – Fundamentals and* 

[40] A. M. Kern and O. J. F. Martin, "Surface integral formulation for 3D simulations of plasmonic and high permittivity nanostructures", J Opt. Soc. Am. A, vol. 26, no. 4, pp.

[41] B. Gallinet, A. M. Kern, and O. J. F. Martin, "Accurate and versatile modeling of electromagnetic scattering on periodic nanostructures with a surface integral

approach", J Opt. Soc. Am. A, vol. 27, no. 10, pp. 2261-2271, 2010.

medium, *IEEE Trans. Antennas Propagat.*, vol. 46, pp. 397-406, Mar. 1998.

doi:10.1029/2010RS004582, 21 May 2011.

*Waves*, vol. 25, no. 8. pp. 1263-1270, 2004.

manual, 2006, http://www.cst.com

vol. 46, pp. 407-413, March 1998.

*Techniques*, vol. 50, pp. 155-164, January 2002.

*Waves*, vol. 25, no. 8. pp. 1263-1270, 2004.

*Applications*, vol. 8, no. 4, pp. 278-284, 2009.

*Techniques*, vol. 52, no. 1, pp. 175-182, January 2004.

[27] www.ansoft.com [28] www.comsol.com [29] www.jcmwave.com

magmas/

85-94, 2010.


46 Plasmonics – Principles and Applications

48, Dec. 1996.

Artech House, 1997.

1086-1094, April 2007.

[15] www.lumerical.com

1966.

1968.

4, Issue 18, pp.387 – 389, Sept. 1968.

Press., 2000, ISBN 978-0-7803-4747-2.

Electromagnetics., pp. 30-34, Nov. 1991.

Journal, Vol. 4, No. 1, pp. 267-282 , Jan. 2012.

*and Theoretical Nanoscience*, vol. 6, no. 3, pp. 763-774, 2009.

[8] J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite element method for

[9] R. L. Courant, "Variational methods for the solution of problems of equilibrium and

[10] S. Ahmed, "Finite-element method for waveguide problems," Electronics Letters, vol.

[11] R. Coccioli, T. Itoh, G. Pelosi, P. P. Silvester, "Finite-Element Methods in Microwaves: A Selected Bibliography," IEEE Antennas and Propag. Magazine, vol. 38, Issue 6, pp.34 –

[12] K. S. Yee, "Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media," IEEE Trans on Antennas Propag, vol 14, pp. 302-307,

[13] Taflove, Computational electrodynamics: the finite difference time domain method,

[14] D. M. Sullivan, "Electromagnetic simulation using the FDTD method", Wiley – IEEE

[16] T. Weiland, "A discretization method for the solution of Maxwell's equations for sixcomponent fields", Electronics and Communications AEÜ, vol 31, No. 3, 116–120, 1977. [17] Bossavit and L. Kettunen,"Yee-like schemes on a tetrahedral mesh, with diagonal

[18] M. Celuch-Marcysiak and W. K. Gwarek, "Comparative study of the time-domain methods for the computer aided analysis of microwave circuits," Int. Conf. on Comp. in

[19] R. F. Harrington, "Field Computation by Moment Methods", New York: Macmillan,

[20] T. K. Sarkar, E. Arvan, and S. Ponnapalli, "Electromagnetic scattering from dielectric

[21] Y. Schols and G. A. E. Vandenbosch, "Separation of horizontal and vertical dependencies in a surface/volume integral equation approach to model quasi 3-D structures in multilayered media ", IEEE Trans. Antennas Propagat., vol. 55, no. 4, pp.

[22] M. Vrancken and G.A.E. Vandenbosch, Hybrid dyadic-mixed-potential and combined spectral-space domain integral-equation analysis of quasi-3-D structures in stratified

[23] X. Zheng, V. K. Valev, N. Verellen, Y. Jeyaram, A. V. Silhanek, V. Metlushko, M. Ameloot, G. A. E. Vandenbosch, and V. V. Moshchalkov, "Volumetric Method of Moments and conceptual multi-level building blocks for nanotopologies", Photonics

[24] J. Smajic, C. Hafner, L. Raguin, K. Tavzarashvili, and M. Mishrickey, "Comparison of Numerical Methods for the Analysis of Plasmonic Structures", *Journal of Computational* 

bodies", IEEE Trans. Antennas Propag., vol. 37, pp. 673-676, May 1989.

media, IEEE Trans Microwave Theory Tech, vol 51, pp. 216-225, 2003.

vibration," Bulletin of the American Mathematical Society, 5, pp. 1-23, 1943.

electromagnetics, IEEE Press, Oxford University Press, 1997.

lumping," Int. J. Numer. Model., vol. 42, pp. 129 – 142, 1999.


[42] I. Chremmos, "Magnetic field integral equation analysis of surface plasmon scattering by rectangular dielectric channel discontinuities", J Opt. Soc. Am. A Opt. Image Sci. Vis., vol. 27, no. 1, pp. 85-94, 2010.

**Chapter 3** 

© 2012 Gao, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

© 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

distribution, and reproduction in any medium, provided the original work is properly cited.

and reproduction in any medium, provided the original work is properly cited.

**Numerical Simulations** 

Xingyu Gao

http://dx.doi.org/10.5772/50751

**1. Introduction** 

technologies.

**of Surface Plasmons Super-Resolution** 

Surface plasmons(SPs) are evanescent waves that propagate along the surface of dispersive media [1]. When a light beam incidents onto the surface of metal, the electromagnetic field interacts with the free electrons inside the metal, which leads to the oscillation of free electrons to excite the SPs in the exit medium. Noble metals with nanostructure geometry have special optical properties because they can excite localized surface plasmons (LSPs) under the illumination of light field. Since the pioneer study of Ritchie [2], the special optical properties of SPs in nano-scale have been tightly investigated. Various metallic nanostructures are reported for nano-plasmonic devices, including thin film[3,4], nanowires[5,6], nanorods[7,8], nano-hole array[9,10], nano-slits[11] et. al. Due to the sub-wavelength excitation range and strong field enhancement properties of SPs, they have been widely applied in super-resolution optical microscopy[12], nano-photonic trapping technology[13],

In this chapter we firstly introduce a numerical algorithm for implementing the relative permittivity model of dispersive media in finite difference time domain(FDTD) method, i.e. piecewise linear recursive revolution(PLRC) method. Next the super-resolution phenomenon derived from the surface plasmons excited at the focal region is presented and analyzed. In the third part, we explore two kinds of plasmonic nano-waveguide: parallel nanorods and metal-dielectric-metal structure. Their optical properties, such as long distance waveguiding from the focal region, turning waveguiding effect and optical switch effect, will be demonstrated in detail. Finally, we conclude our research results in recent years and look forward the applications of surface plasmons in nano-

**Focusing and Nano-Waveguiding** 

biology and medical sciences[14] and nano-photonic waveguide[15].

Additional information is available at the end of the chapter

