合成在特定波長有表面電漿共振之金屬奈米結構

(1)

國立臺灣大學電機資訊學院光電工程學研究所碩士論文

Graduate Institute of Photonics and Optoelectronics College of Electrical Engineering and Computer Science

National Taiwan University Master Thesis

合成在特定波長有表面電漿共振之金屬奈米結構 Synthesis of Metallic Nanostructures with Surface-Plasmon

Resonance at Designated Wavelengths

歐陽良昱 Liang-Yu Ou Yang

指導教授：江衍偉博士

中華民國98年7月 July 2009

(2)

(3)

誌謝

『桃李春風一杯酒，江湖夜雨十年燈』，碩士兩年，首先感謝江

衍偉老師的耐心指導，不論是在研究方面或電磁學，都讓我收穫良

多。老師對於學問的純真追求，堪稱典範，跟江老師討論電磁問題，

是我在碩士生涯中最快樂的時光。另外也感謝張宏鈞、吳育任、楊志

忠老師對論文提出寶貴意見。接著要感謝我的父母，他們是我人生中

最重要的人，有他們在經濟與心靈上的支持，讓我有勇氣去面對任何

難題。也感謝郭仰、王志洋、蔡富吉、莊文宏、林俊安、陳弘履、沈

昌煒等學長在理論與程式上的幫助，讓我的研究能順利完成。方嘉

鵬、吳碩彥同學則是在學習、玩樂、修課、聯誼上不可多得的好夥伴，

而于贊堯、廖家鴻、張然凱學弟的加入，對實驗室又注入了一份熱情

與活力。最後，要感謝禎禧，在這兩年裡默默的陪伴著我。

歐陽良昱謹識

予台大光電工程研究所

民國九十八年七月

(4)

摘要

我們利用模擬退火法與邊界積分方程法來合成在特定波長具有

表面電漿共振特性的金屬奈米結構。數值結果依應用可分為: 藍光

LED (波長 435 nm)、白光 LED (波長 450 nm 與 570 nm)、與綠光 LED

(波長 535 nm)。藍光 LED 係由二金屬圓柱部分嵌於金屬半空間之奈

米結構所組成。吾人首先合成一結構，其表面電漿極化子和局域表面

電漿子在單一波長 (435 nm) 能有效耦合。接著再合成另一結構，使

其在另一波長 (520 nm) 僅有局域表面電漿子之共振。數值模擬顯示:

表面電漿極化子和局域表面電漿子在單一波長有效耦合時，確實可提

升其附近電偶極之輻射量。白光 LED 係由二分離之金屬圓柱與一金屬

半空間之奈米結構所組成。吾人合成一結構，在藍光波長 (450 nm)

與黃光波長 (570 nm) 皆具有局域表面電漿子共振，因而提升混光而

成的白光發光量。綠光 LED 係由二分離之金屬橢圓柱與一金屬半空間

之奈米結構所組成。吾人合成一結構，其表面電漿極化子和局域表面

電漿子在單一波長 (535 nm) 能有效耦合，因而提升電偶極之輻射量

與輻射效率。

(5)

Abstract

The simulated annealing and the boundary integral-equation

methods are used to synthesize metallic nanostructures with

surface-plasmon resonance properties at designated wavelengths. The

numerical results include three parts according to the possible

applications: blue-light LEDs (wavelength 435 nm), white-light LEDs

(wavelengths 450 nm and 570 nm), and green-light LEDs (wavelength

535 nm). For the blue-light LEDs, the structure consists of two metallic

circular cylinders partially embedded in a metallic half-space. We first

synthesize a metallic nanostructure such that the surface plasmon

polariton (SPP) and the localized surface plasmon (LSP) couple

effectively at their common resonant wavelength (435 nm). Next, we

synthesize another structure for optimization at wavelength 520 nm, at

which only the LSP resonance occurs. From numerical simulations, it is

demonstrated that the enhancement of the dipole emission is better for

optimization at wavelength 435 nm than that at wavelength 520 nm. In

the aspect of white-light LEDs, the structure is composed of two separate

metallic circular cylinders and a metallic half-space. We synthesize a

(6)

metallic nanostructure, which has LSP resonances at wavelength 450 nm

(blue light) and at wavelength 570 nm (yellow light), leading to the

enhancement of white-light emission. For the green-light LEDs, the

structure consists of two separate metallic elliptical cylinders and a

metallic half-space. We synthesize a metallic nanostructure such that the

SPP and LSP couple effectively at their common resonant wavelength

(535 nm), leading to the enhancement of both the dipole emission and the

emission efficiency.

(7)

Chapter 1 Introduction

...1

Chapter 2 Surface Plasmon (SP)

...4

2.1 Surface plasmon polariton (SPP)...4

2.2 Localized surface plasmon (LSP)...6

2.3 Dipole-SP coupling phenomenon...8

Chapter 3 Numerical Methods

...13

3.1 Boundary integral-equation method (BIEM) ...13

3.2 Simulated annealing method ...18

3.2.1 Basic principles...18

3.2.2 Synthesis of metallic nanostructures by using simulated annealing method...20

Chapter 4 Numerical Results

...24

4.1 Metallic nanostructures consisting of two circular cylinders partially embedded in a half-space ...24

4.2 Metallic nanostructures consisting of a half-space and separate circular cylinders...31

4.3 Metallic nanostructures consisting of a half-space and separate elliptical cylinders ...39

Chapter 5 Conclusions

...91

References

...93

(8)

List of Figures

2.1 Electric field distributions of a SPP at a metal-dielectric flat interface ...11 2.2 SPP dispersion relation at a metal-dielectric flat interface ...11 2.3 Schematic diagram of the electron-hole recombination and QW-SP

coupling ...12 3.1 Schematic diagram of a single-object scattering problem treated by

the BIEM ...22 3.2 (a) A boundary for applying the equivalent electric and magnetic

surface currents in a 2D TM case problem. (b) The discretized

version of (a)...22 3.3 Flow chart of a simulated annealing process ...23 4.1 Variation of downward emission during the iteration process...55 4.2 Spectra of downward emission and enhancement factor of

downward emission for three different structures...56 4.3 Distributions of the absolute values of magnetic field for

structure A at wavelength 435 nm (a), 404 nm (b), 465 nm (c), and

520 nm (d) ...57 4.4 Distributions of the absolute values of electric field for structure A

at wavelength 435 nm (a), 404 nm (b), 465 nm (c), and 520 nm (d) ...58 4.5 Distributions of the absolute values of magnetic field (a) and

electric field (b) for structure B at wavelength 520 nm ...59 4.6 Distributions of the absolute values of magnetic field (a) and

electric field (b) for structure D at wavelength 435 nm ...59 4.7 Variation of cost function in the iteration process ...60 4.8 Variations of downward emission at two wavelengths 450 nm and

570 nm in the iteration process...61 4.9 Spectra of total emission (dashed curve) and downward emission

(solid curve) for structure E...62 4.10 Spectra of enhancement factor of total emission (dashed curve)

and downward emission (solid curve) for structure E. ...62 4.11 Distributions of the absolute values of magnetic field for structure

E at wavelength 435 nm (a), 450 nm (b), 540 nm (c), and 570 nm

(d)...63

(9)

4.12 Distributions of the absolute values of electric field for structure E

at wavelength 435 nm (a), 450 nm (b), 540 nm (c), and 570 nm (d) ...64

4.13 Spectra of total emission (dashed curve) and downward emission (solid curve) for structure F...65

4.14 Spectra of enhancement factor of total emission (dashed curve) and downward emission (solid curve) for structure F...65

4.15 Distributions of the absolute values of magnetic field for structure F at wavelength 435 nm (a), 450 nm (b), 470 nm (c), 540 nm (d), and 570 nm (e)...66

4.16 Distributions of the absolute values of electric field for structure F at wavelength 435 nm (a), 450 nm (b), 470 nm (c), 540 nm (d), and 570 nm (e)...67

4.17 Spectra of total emission (dashed curve) and downward emission (solid curve) for structure G ...68

4.18 Spectra of enhancement factor of total emission (dashed curve) and downward emission (solid curve) for structure G ...68

4.19 Distributions of the absolute values of magnetic field for structure G at wavelength 435 nm (a), 450 nm (b), 540 nm (c), and 570 nm (d)...69

4.20 Distributions of the absolute values of electric field for structure G at wavelength 435 nm (a), 450 nm (b), 540 nm (c), and 570 nm (d) ...70

4.21 Variation of cost function in the iteration process ...71

4.22 Variation of downward emission in the iteration process ...72

4.23 Variation of emission efficiency in the iteration process...72

4.24 Spectra of total emission (dashed curve) and downward emission (solid curve) for structure H ...73

4.25 Spectra of enhancement factor of total emission (dashed curve) and downward emission (solid curve) for structure H ...73

4.26 Spectrum of emission efficiency for structure H ...74

4.27 Distributions of the absolute values of magnetic field for structure H at wavelength 435 nm (a) and 535 nm (b)...74

4.28 Distributions of the absolute values of electric field for structure H at wavelength 435 nm (a) and 535 nm (b)...75

4.29 Spectra of total emission (dashed curve) and downward emission (solid curve) for structure I...75

4.30 Distributions of the absolute values of magnetic field for structure I at wavelength 535 nm (a) and 555 nm (b) ...76 4.31 Distributions of the absolute values of electric field for structure I

(10)

at wavelength 535 nm (a) and 555 nm (b)...76

4.33 Variation of downward emission in the iteration process ...78

4.34 Variation of emission efficiency in the iteration process...78

4.35 Spectra of total emission (dashed curve) and downward emission (solid curve) for structure J...79

4.36 Spectra of enhancement factor of total emission (dashed curve) and downward emission (solid curve) for structure J ...79

4.37 Spectrum of emission efficiency for structure J...80

4.38 Distributions of the absolute values of magnetic field for structure J at wavelength 415 nm (a), 435 nm (b), and 440 nm (c) ...81

4.39 Distributions of the absolute values of electric field for structure J at wavelength 415 nm (a), 435 nm (b), and 440 nm (c)...82

4.41 Spectra of total emission (dashed curve) and downward emission (solid curve) for structure K. ...84

4.42 Spectra of enhancement factor of total emission (dashed curve) and downward emission (solid curve) for structure K ...84

4.43 Spectrum of emission efficiency for structure K ...85

4.44 Distributions of the absolute values of magnetic field (a) and electric field (b) for structure K at wavelength 435 nm ...85

4.46 Spectra of total emission (dashed curve) and downward emission (solid curve) for structure L...87

4.47 Spectrum of emission efficiency for structure L ...87

4.48 Distributions of the absolute values of magnetic field for structure L at wavelength 435 nm (a) and 535 nm (b) ...88

4.49 Distributions of the absolute values of electric field for structure L at wavelength 435 nm (a) and 535 nm (b)...88

4.50 Geometry of structure M ...89

4.51 Variations of total emission (dashed curve) and downward emission (solid curve) for structure M at wavelength 535 nm with the horizontal position of the dipole source ...89

4.52 Variations of enhancement factor of total emission (dashed curve) and downward emission (solid curve) for structure M at wavelength 535 nm with the horizontal position of the dipole source...90

(11)

List of Tables

4.1 Peak enhancement factors of downward emission for three structures in a spectral range from 400 nm to 700 nm ...54 4.2 Parameters for all structures discussed in section 4.3 (except for

structure M) and their downward emission peak values and the

corresponding wavelengths ...54

(12)

Chapter 1 Introduction

The phenomenon of surface plasmon (SP) has been explored since the beginning of last century. The SP represents the behavior which electrons collectively oscillate near the surface of metal. For the applications of SP in the optical field and nanophotonics, the devices using SP in sub-wavelength scale enable us to miniaturize integrated optical circuits for confining or controlling electromagnetic wave. SP waveguides are typical examples. In chemical areas, the SP can be applied to enhance the intensity of signal for Raman scattering. This technique is called “surface enhanced Raman spectroscopy (SERS)”. In biological areas, the idea of SP field enhancement is used to construct the biosensors for detecting interesting molecules with high sensitivity. In recent years, the coupling of SP with a light emitter for enhancing the emission efficiency has been widely studied [1-18]. The reason for emission enhancement can be explained by the classical electrodynamics.

If the SP mode has a large field value in the neighborhood of a dipole, the radiation of that dipole will be enhanced [19]. It means the enhancement

(13)

of spontaneous emission rate through the coupling between quantum wells and SP is possible. Hence, the SP-assisted devices to improve light extraction efficiency and internal quantum efficiency have become attractive. Many approaches have been proposed to enhance the light extraction efficiency, such as LEDs with resonant cavity, surface texturing, and output coupling surface plasmon at corrugated or grating metallic surfaces.

In general, the SP modes can be classified as the surface plasmon polariton (SPP) and localized surface plasmon (LSP). In this thesis, we focus on discussing the coupling properties of both SPP and LSP. The coupling of SPP and LSP are numerically investigated in the frequency domain. The organization of this thesis is outlined as follows.

In chapter 2, the fundamentals of the SP are briefly reviewed. In chapter 3, we introduce our simulation methods, including boundary integral-equation method (BIEM) and simulated annealing (SA) technique. Compared to the commonly used finite-difference time-domain (FDTD) method, the BIEM can deal with complex geometries easier than FDTD. The reason is that FDTD is a formulation using fixed square space grids, but BIEM is with unstructured-mesh

(14)

based on a surface formulation, which automatically takes care of discontinuities between materials and models arbitrarily shaped structures.

The SA is a probabilistic and iterative method of optimization mimicking the real annealing process of matter [20, 21]. We also introduce the theory of SA and its application to synthesizing metallic nanostructures in this chapter. In chapter 4, the numerical simulation results are presented.

Finally, some conclusions are drawn in chapter 5.

(15)

Chapter 2 Surface Plasmon (SP)

2.1 Surface plasmon polariton (SPP)

A plasma is a quasi-neutral gas of the electrons, positive ions, and neutral particles possessing certain collective behaviors. The tendency to restore the charge neutrality leads to the plasma oscillation with an angular frequency ω_p given by

2

0 p

e

Ne ω m

= ε , (2.1) where N is the number density of electrons, m_e is the mass of an electron, ε₀ is the permittivity in free space, and e is the electric charge of an electron.

A noble metal, such as Ag or Au, can be regarded as a cold plasma.

In other words, we can consider it as a dielectric with an effective dielectric constant given by the Drude model [22]

) 1 (

) (

2

γ ω ω ω ω

ε i

p

m = − + , (2.2) where γ is the damping constant or collision frequency. In this study, the dielectric constant of Ag is assumed to follow the Drude model with the angular plasma frequency at ωp = 1.19×10¹⁶(rad·s^-1) and the damping

(16)

constant at γ = 1.32×10¹⁴(rad·s^-1). For another metal Au, ωp = 9.37×10¹⁵ (rad·s^-1) and γ = 1.04×10¹⁴(rad·s^-1).

In the bulk of a cold plasma, the plasma oscillation can be described as an oscillatory fluctuation of the electron density. On the other aspect, there is another mode distinct from the bulk plasma mode, known as

“surface plasmon”. It can exist on the interface between a metal and a dielectric, and is associated with the oscillation of surface charges excited by exterior electromagnetic field. Note that the electromagnetic field and polarization wave can be quantized with energy quanta called polaritons.

Hence, a combination of excitation of surface plasmon and polariton is called the surface plasmon polariton (SPP). A SPP is evanescent perpendicular to the metal-dielectric interface. This nonradiative surface wave can be analyzed by solving Maxwell’s equations and matching the associated boundary conditions.

Figure 2.1 shows the interaction between the electromagnetic wave and the surface plasma oscillation at a metal-dielectric flat interface. In the figure, we consider the p-polarized (or TM) wave with the electric field lying in the plane of incidence. It can be shown that the x-component wavevector of SPP is given by [23]

(17)

( ) ( )

m d

SPP

m d

k c

ε ω ε ω

ε ω ε

= + . (2.3)

Here, ε ω_m( ) is the frequency-dependent dielectric constant of the metal and ε_d is the dielectric constant of the dielectric.

For the lossless Drude model, k_SPP approaches infinity as ω tends to the SPP resonant frequency

d p

SPP ε

ω ω

= +

1 . (2.4) Figure 2.2 shows the dispersion relation of SPP at a metal-dielectric flat interface. Here, the solid line represents the SPP dispersion curve, and the dashed one shows the dispersion curve of the propagating wave in the dielectric medium, known as the light line. Note that the dispersion curve of a SPP mode is located to the right of the dielectric light line, showing the momentum mismatching. This is the problem that a SPP cannot be excited by the direct illumination of light from the dielectric side. Besides, as ω approaches ω_SPP, the group velocity of SPP mode goes to zero.

Since the decay rate is very high in either dielectric or metal region, a SPP form localized fluctuation at a metal-dielectric interface.

2.2 Localized surface plasmon (LSP)

The localized surface plasmon (LSP) occurs near the surface of some

(18)

metallic nanoparticles. The resonance of LSP depends mainly on the structure geometry and the polarization state of the incident light.

Especially for nanoscale metallic geometry, such as a metallic sphere or a metallic circular cylinder, the LSP excitation is accompanied with a highly localized field around the nanostructure. Let us consider a metallic circular nano-cylinder in a homogeneous, linear, isotropic dielectric medium, under the illumination of visible light. Since the wavelength is much larger than the dimension of the circular cylinder, the quasi-static approximation is applicable to this problem. By solving Maxwell’s equations and matching the associated boundary conditions, we can derive the formula for the LSP resonant frequency [24]

d cylinder p

circular

LSP ε

ω ω

= +

1 . (2.5) For a single Ag circular cylinder embedded in GaN, the LSP resonant wavelength is at 430 nm. Note that the dielectric constant of GaN is set to be 6.25 in all the following simulations.

When the quasi-static approximation is not appropriate, the full Maxwell’s equations have to be solved. For a sphere, the Mie theory has been well known [25]. For a circular cylinder, similar discussion can be

(19)

Note that, in the microwave region, the metal can almost be regarded as a perfect electric conductor (PEC), and the phenomenon of SP is unclear. The SP mainly occurs in the visible spectral range for noble metals.

2.3 Dipole-SP coupling phenomenon

Okamoto et al. proposed possible mechanisms of quantum well (QW)-SP coupling and light extraction, as shown in Fig. 2.3 [1-3]. First, excitons are generated in the QW by photo-pumping or electrical-pumping. For a sample uncoated with metal, these excitons are terminated by the radiative (P_rad) or nonradiative (P_non) recombination rates, and the internal quantum efficiency (η_int) is defined as

int P_rad /(P_rad P_non)

η = + . (2.6) When a metallic layer is grown within the near-field of the active layer, and when the bandgap energy (ω_BG) of InGaN active layer is close to the electron vibration energy (ω_SP) of SP at the metal-semiconductor surface, then the QW energy can transfer to the SP. The photoluminescence (PL) decay rate is enhanced through the QW-SP coupling rate (P_SP), which is expected to be very fast. High electromagnetic field is produced by the large density of states from the SP dispersion diagram, making P_SP

(20)

increase. The QW-SP coupling in LED devices may be considered beneficial to the light extraction efficiency, because the SP can re-couple to radiate light by some specified fabrication. If the metallic surface is perfectly flat, the SP energy would be thermally dissipated without radiation. However, the SP energy can be extracted as light by providing roughness or nano-structuring the metallic layer. Such roughness provides SP additional momentum to scatter, to lose momentum, and to couple to radiation. The few tens of nanometer sized roughness in the Ag surface layer can be obtained by controlling the evaporation conditions or by microfabrication to obtain the high photon extraction efficiency.

In order to understand the strength of QW-SP coupling, Fermi’s golden rule has been applied to estimate the transition probability or the recombination rate, Γ_p( )ω as

1 2 2

( ) ( ) ( )

p ω π f d E r ie ρ ω Γ = =τ

i

. (2.7) Here, τ is the spontaneous decay constant, f d E r i( )_e

i is the dipole interaction matrix element and ρ ω( ) is the plasmon mode density of states (DOS) [27, 28]. The density of states ρ ω( ) can be obtained from the derivative

( 2) d k

dω of the dispersion relation ω( )k for a SP mode.

Then Eq. (2.7) becomes

(21)

2 2

2 2 ( )

( ) ( )

4 ( )

p e

L d k f d E r i

d ω π

π ω

 

Γ =  

 

i

. (2.8) By increasing the density of states of a SP mode, the spontaneous recombination rate can be enhanced as the photon frequency of the active layer approaches the SP resonance frequency. Besides, when the active layer is closer to the metal-dielectric surface, the QW-SP coupling becomes stronger.

(22)

Fig. 2.1 Electric field distributions of a SPP at a metal-dielectric flat interface.

Fig. 2.2 SPP dispersion relation at a metal-dielectric flat interface.

(23)

Fig. 2.3 Schematic diagram of the electron-hole recombination and QW-SP coupling.

(24)

Chapter 3 Numerical Methods

3.1 Boundary integral-equation method (BIEM)

The boundary integral-equation method used here is a 2D version of the surface integral-equation method based on the Stratton-Chu formulas [29-35]. The Stratton-Chu formulas state that the total field at an observation point ρ, in an enclosed homogeneous, (^{ε µ}^, ), region can be calculated as the sum of incident field and scattered field. The scattered field is produced by the equivalent electric and magnetic surface currents,

(^{J M}^, ), on all boundaries of this region. The expressions of the total fields are given by

( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

+ , , ,

inc

E E

i J M i J dl

H H

i M J i M dl

ρ ρ

ωµ ρ ϕ ρ ρ ρ ϕ ρ ρ ρ ϕ ρ ρ

ωε

ρ ρ

ωε ρ ϕ ρ ρ ρ ϕ ρ ρ ρ ϕ ρ ρ

ωµ

 =

  ′ ′ − ′ ×∇′ ′ − ∇ ⋅′ ′ ∇′ ′  ′

    



 =

    

  ′ ′ + ′ ×∇′ ′ − ∇ ⋅′ ′ ∇′ ′  ′

    



∫

.

(3.1) Here, ϕ ρ ρ( , ′)

is the 2D Green’s function of the homogeneous region given by

(25)

( , ) 0^{( )}¹

( )

4 i H k

ϕ ρ ρ ′ = ρ ρ− ′

, (3.2)

where H₀^{( )}¹ is the first-kind Hankel function of order zero, ω is the angular frequency, and k =ω µε . Note that the time-harmonic variation of form e⁻^{i t}^ω is assumed in this thesis.

If the observation point is located on a boundary (approaching from the interior of the region), then the electric and magnetic fields become

( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

2

+2 , , ,

2

+2 , , ,

inc

E E

i J M i J dl

H H

i M J i M dl

ρ ρ

ωµ ρ ϕ ρ ρ ρ ϕ ρ ρ ρ ϕ ρ ρ

ωε

ρ ρ

ωε ρ ϕ ρ ρ ρ ϕ ρ ρ ρ ϕ ρ ρ

ωµ

 =

  ′ ′ − ′ ×∇′ ′ − ∇ ⋅′ ′ ∇′ ′  ′

    



 =

    

  ′ ′ + ′ ×∇′ ′ − ∇ ⋅′ ′ ∇′ ′  ′

    



∫

.

(3.3) Note that the above line integrals are evaluated as the Cauchy principal values. Once the field expressions on each side of an interface approaching from different regions expressed by Eq. (3.3) are derived, we can match their tangential components to obtain a set of boundary integral equations associated with that interface.

Figure 3.1 shows a schematic diagram of a typical problem with a simple geometry containing only two regions. Region 2 is an infinitely extended background containing current sources (^J^inc^,^M^inc) ^which

(26)

produce the incident fields, while region 1 is a scatterer without any source inside it. The boundary integral equations on interface C

between the two regions can be derived as [35]

( ) ( ) ( )

( ) ( )

2 1

2 2 1 1 2 1

2 1

2 2 1 1 2 1

2 1

ˆ ˆ

inc

s

inc

s

n E

n i J M i J dl

n H

n i M J i M dl

φ φ

ω µ φ µ φ φ φ

ω ε ε

φ φ

ω ε φ ε φ φ φ

ω µ µ

 ×

  

 = × − + + × ∇′ + + ∇ ⋅′ ∇′ +  ′

  

  

 ×

   

 = × − + − × ∇′ + + ∇ ⋅′ ∇′ +  ′

   



∫

,

(3.4) where we define n^ˆ as the unit normal vector of the boundary C

pointing to R₂. Note that Eq. (3.4) is obtained by matching the tangential electric and magnetic fields on the two sides of C, with the fields expressed by Eq. (3.3).

To numerically solve the integral equation, one has to discretize the original problem into a matrix form. A commonly used technique is the method of moments (MoM), which expands the unknown equivalent surface currents (^{J M}^, ) with a set of local bases. The basis-expanded integral equation is further tested with the local bases (by Galerkin testing) and then transformed into a matrix equation.

In our formulation, the linear basis is used rather than the pulse basis,

(27)

problems. Figure 3.2 (a) shows a boundary, on either side of which the equivalent surface currents, (^{J M}^, )^and (⁻^J^,⁻^M ), are to be expanded with the linear bases. Here, a 2D transverse magnetic (TM) polarization, defined with field components E_x, E_y and H_z, is assumed. Figure 3.2 (b) is a discretized version of Fig. 3.2 (a). In Fig. 3.2 (b), the original smooth curve is approximated by several piecewise lines, with the nodes

ρn used to define the linear basis. Two sets of linear bases f_j

and g_j

are used for expanding J

and M

, respectively. Both f_j

and g_j

consist of two parts, i.e., f_j =t^ˆ_j_,1f_j_,1+t^ˆ_j_,2f_j_,2

and ^g^j ⁼^{z g}^ˆ

(

^j^,1⁺^g^j^{, 2}

)

^{, where}

ˆ,1

tj and t^ˆ_j_,2 are unit tangent vectors defined by

( )

,1 1 1

,2 1 1

ˆ

j j j j j

t

ρ ρ ρ ρ

− −

+ +

 = − −



 = − −



, (3.5)

so that t^ˆ_j_{,1 2}_{( )} is always pointed to the direction of ascending label number of the nodes. Also given below are the definitions of f_j_,1 and f_j_,2

1

,1 1

,

0, otherwise

j

j j

j j j

f

ρ ρ

ρ ρ ρ ρ ρ

−

 ′ −

′∈

= −



, ¹ ¹

,2 1

,

0, otherwise

j

j j

j j j

f

ρ ρ

ρ ρ ρ

ρ ρ

+

+ +

 − ′

′∈

= −



. (3.6)

Note that g_j_{,1 2}_{( )} = f_j_{,1 2}_{( )} is also defined.

Here a convention of the equivalent surface currents should be stated.

We define the equivalent surface currents with positive sign in the following manner

(28)

( )

,1 ,1 ,2 ,2

,1 , 2

ˆ ˆ

ˆ

j j j j j

j

j j j

j

J J t f t f

M z M g g

 = +



= +



∑

. (3.7)

Those with negative sign are denoted by (⁻^J^,⁻^M) , with their definitions given below

( )

,1 ,1 ,2 ,2

,1 ,2

ˆ ˆ

ˆ

j j j j j

j

j j j

j

J J t f t f

M z M g g

− = − +



− = − +



∑

. (3.8)

Note that the directions of the equivalent surface currents plotted in Fig.

3.2 (a) do not really mean the true directions of the equivalent surface currents. The true directions of them are determined by the unknown coefficients, J_j and M_j, and the fixed unit vectors, t^ˆ_j_,1, t^ˆ_j_,2, and ˆz. Then the Galerkin testing procedure is used to execute the field expression on the interface, and we can transform the boundary integral equation into a matrix equation with the form of

[ ]

^A^x⁼^b, Here the all elements of matrix

[ ]

^A are known after testing, column vector b

represents the incident field, and column vector x

represents unknown electric and magnetic currents. Once the equivalent surface currents are obtained through matrix inversion, the complex electromagnetic field at any position can be readily calculated through the boundary integral.

(29)

3.2 Simulated annealing method

How to find an optimal solution is an old problem, and it exists in various areas. Many mathematical techniques have been proposed to treat optimization problems [20, 21]. The simulated annealing (SA) method, which is related to the Monte Carlo method, is a useful algorithm for treating optimization problems. It not only has good efficiency for local searching but also is capable of searching a global extremum among several possible local extrema. The basic principles of the SA and how we apply it to synthesize metallic nanostructures will be described in the following.

3.2.1 Basic principles

The theory of SA algorithm stemmed from the pioneering work by Metropolis in 1953 [21]. In 1983, Kirkpatrick proposed this algorithm for solving the combinatorial optimization problem [20]. The SA is a probabilistic and iterative method of optimization mimicking the real annealing process of matter from the liquid state to the solid state.

Initially, the temperature is high and the atoms of the liquid material are unstable in such high energy states. When the temperature is slowly decreased, the atoms will be gradually stabilized. If the temperature is

(30)

low enough, the material will become a perfect crystal in this low energy state. Similarly, an optimization problem can be posed as to find a state at which the cost function (analogous to the energy) attains a minimum.

To start with SA, we need to choose certain operation quantities, including an initial state and a cost function C. The latter changes its value at different states and is to be minimized. Also an appropriate temperature function T, whose value decreases with time, is selected.

During the optimization process, the state is subject to a random change such that a new value of the cost function C is obtained. The probability, P, of choosing the new state to replace the old one is determined by ∆C, the difference between the new and old values of the cost function, as

{

1 , if 0 , if 0

C T

C

e C

P

= _−∆ ∆ <

∆ ≥ ^, (3.9) where the exponential form is adopted for initiating the Boltzmann distribution in thermodynamics. Based on Eq. (3.9), a better state of a lower cost value, i.e., ∆C <0, is accordingly preferred. When a worse state with the condition deviating away from the target occurs, i.e.,

≥0

∆C , the selection is still implemented if ^exp(⁻^{C T}^{/ )} is larger than a computer-generated random number between 0 and 1. In other words,

(31)

there is certain probability that we accept a worse state to avoid being trapped in a local minimum. Therefore, we reject the change of state if the cost function rises and the above judgment is not satisfied. After the evolution process mentioned above, the temperature T is decreased. We iterate such steps until the designated criterion is met. Figure 3.3 shows the flow chart of a general SA process.

Note that the efficiency of SA seriously depends on the initial state and the temperature profile. If the initial state is quite different from the optimal state, the iteration would be time-consuming. If the temperature decreases too rapidly, the worse result would hardly be accepted during the iteration and the state may be trapped in a local minimum. However, the slower the temperature decreases, the more iteration steps will be needed to achieve the optimization. Therefore, properly controlling the temperature and selecting the initial state would be crucial for obtaining a satisfactory optimization result efficiently.

3.2.2 Synthesis of metallic nanostructures by using simulated annealing method

Now, the SA is used to synthesize metallic nanostructures. In our simulation, the adjustable parameters include the radius r of each circular

(32)

cylinder (For each elliptical cylinder, we choose the length a of major semiaxis to adjust, and the length of minor semiaxis is fixed), the separation d between centers of the two circular or elliptical cylinders, the distance t from the center of each circular or elliptical cylinder to the plane of flat metal, and the vertical displacement h from the plane of flat metal to the dipole source. But in section 4.1, parameter h is fixed, unless specified otherwise. In every iteration, the probability of randomly adjusting these parameters by + 2 nm, - 2 nm, or 0 nm is equal. However, we randomly change parameter t by 0, ±

10

r with equal possibility in

section 4.1. The cost function, in general, consists of two factors and the proportion of these two factors will be adjusted to keep their balanced weights in the cost function. Nevertheless, the cost function in section 4.1 is only as the negative value of downward emission.

(33)

Fig. 3.1 Schematic diagram of a single-object scattering problem treated by the BIEM.

(a) (b)

Fig. 3.2 (a) A boundary for applying the equivalent electric and magnetic surface currents in a 2D TM case problem. (b) The discretized version of (a).

(34)

Fig. 3.3 Flow chart of a simulated annealing process.

(35)

Chapter 4 Numerical Results

4.1 Metallic nanostructures consisting of two circular cylinders partially embedded in a half-space

As mentioned in chapter 2, the SP modes can be classified as the SPP and LSP. It has been demonstrated that either SPP or LSP can individually enhance the dipole emission at its resonant wavelength [1-3].

However, the possibility of effective coupling between SPP and LSP at a common wavelength has not been discussed yet. Note that under the quasi-static approximation [23, 24], the LSP resonant wavelength of a single metallic circular cylinder is equal to the SPP resonant wavelength of a metallic half-space.

In this section, we synthesize a nanostructure consisting of two metallic circular cylinders partially embedded in a metallic half-space by the simulated annealing (SA) method. The nanostructure is so designed that the SPP and LSP couple effectively at their common resonant wavelength (435 nm), leading to the enhancement of the radiation of a nearby dipole source.

The parameter “downward emission” is defined as the power which

(36)

the dipole radiates in the direction away from the metallic structure in all the following simulations. We adjust the structure by the SA to maximize the downward emission at the resonant wavelength, with the cost function being the minus value of downward emission. Then the spectrum of downward emission for the optimal structure is calculated.

The insert of Fig. 4.1 shows the schematic structure which we try to optimize in this section. The structure and the x-directed dipole are infinitely extended along the z axis to form the two-dimensional problem in all the following simulations. Thus, the field components are Ex, Ey, and Hz only. The metal is Ag, and the surrounding dielectric is GaN in this section. Note that the dielectric in synthesized structures is set to be GaN in all the following simulations. All the adjustable parameters and the way of changing the structure are described in section 3.2.2. The plane of flat metal is set at y = 0 and the position of the dipole is set at (x, y) = (0, -h) in the following simulations. Note that we set h = r + t in this section, unless specified otherwise. We synthesize the optimal structures by adjusting these parameters to maximize the downward emission by using the SA method. After an iteration, if the downward emission of the new structure is larger than the old one, we accept the new structure.

(37)

Otherwise, we accept or reject the new structure according to the SA criterion described in chapter 3. It is well known that the SPP resonant frequency of a metallic half-space is

1

p d

ω ε

+ [23], where ε_d is the dielectric constant of the surrounding dielectric. Based on the Drude model, we calculate the SPP resonant wavelength of the Ag-GaN flat interface to be at 430 nm. In our numerical simulation, however, a dipole cannot supply infinite wave numbers, so we choose a slightly longer wavelength at 435 nm for optimization. In order to manifest the advantage of the effective coupling between SPP and LSP, we also consider another structure for optimization at wavelength 520 nm for comparison. Through adjusting variable parameters, the optimal structure at wavelength 520 nm has mainly the LSP effect. We denote the structure optimized at 435 nm as structure A and that optimized at 520 nm as structure B. After optimization by the SA, the parameters of structure A are: r = 10 nm, d = 68 nm, t = 0.8 nm, h = 10.8 nm, and the parameters of structure B are: r = 14 nm, d = 60 nm, t = 7.6 nm, h = 21.6 nm.

Figure 4.1 shows the variation of downward emission during the iteration process. The green curve represents the chosen temperature distribution used in the SA process. The red curve shows the optimization

(38)

process for downward emission at wavelength 435 nm and the blue curve indicates the optimization process for downward emission at wavelength 520 nm. From these curves, we can find that these SA processes indeed converge to the optimal structures. Note that parameter r of structure B is larger than that of structure A. The reason is that as the resonant wavelength red shifts from 435 nm to 520 nm, the size of system should become bigger. There may be two ways for enlarging the size of system.

One is to increase the separation between the two circular cylinders and the other is to increase the size of circular cylinders. Since the circular cylinders cannot be too far away from the source dipole for effectively exciting the LSP effect, the radius of circular cylinders then becomes larger for optimization at a longer wavelength.

Figure 4.2 shows the spectra of downward emission for structures A, B, and C. Here, the parameters of structure C are the same as those of structure A except that the dipole is moved downward by 10 nm, i.e., h = 20.8 nm. The y axis of solid curves is for the downward emission. The y axis of dashed curves is for the enhancement factor of downward emission. The enhancement factor is defined as the downward emission across the surface of y = -(60 + h) nm divided by the similar downward

(39)

emission when the structured Ag-GaN interface in the insert of Fig. 4.1 is replaced by a plane surface between Ag and GaN at y = 0. The black curves represent the spectra of structure A. The blue curves represent the spectra of structure B. The red curves represent the spectra of structure C.

The peak value of downward emission for structure A occurs at wavelength 435 nm with peak value 0.025, and the peak enhancement factor is 6.12. The peak value of downward emission for structure B occurs at wavelength 520 nm with peak value 0.01, and the peak enhancement factor is 2.71. The peak value of downward emission for structure C occurs at wavelength 435 nm with peak value 0.009, and the peak enhancement factor is 2. The SA processes are successful because the peak of downward emission for each structure occurs indeed at the desired wavelength for optimization. We can find that both the peak value of downward emission and the peak enhancement factor of downward emission for structure A are bigger than those for structures B and C. This is consistent with the fact that at the resonant wavelength 435 nm, the SPP and LSP modes can couple effectively. Besides, structure B supports mainly the LSP mode because its resonant wavelength 520 nm is far away from the SPP resonant wavelength. The SPP and LSP modes for

(40)

structure C are both weak because the dipole is too far from the structure to effectively excite the SPP and LSP modes. Obviously, the peak value of downward emission and the peak enhancement factor of downward emission for structure C are less than those for structures A and B.

To further understand the coupling mechanism between SPP and LSP, we plot the distributions of the absolute values of magnetic field in Figs. 4.3 (a)-(d) for structure A at wavelength 435 nm (a), 404 nm (b), 465 nm (c), and 520 nm (d), and plot the distributions of the absolute values of electric field in Figs. 4.4 (a)-(d) for structure A at wavelength 435 nm (a), 404 nm (b), 465 nm (c), and 520 nm (d). The absolute value of electric field is calculated by the following formula in all the following simulations

2 2 y

x E

E

E = +

. (4.1) The SP modes in Figs. 4.3 (a) and 4.4 (a) include both SPP and LSP, and the coupling effect between them is significant because at wavelength 435 nm both SPP and LSP resonate simultaneously. The field distributions are rather weak in Figs. 4.3 (b)-(d) and 4.4 (b)-(d), indicating that neither SPP nor LSP is effectively excited at these wavelengths for structure A.

(41)

Figure 4.5 shows the distributions of the absolute values of magnetic field (a) and electric field (b) at wavelength 520 nm for structure B. By comparing Fig. 4.5 (a) and Fig. 4.3 (d), we find that the LSP effect for structure B is more remarkable than that for structure A at wavelength 520 nm, leading to larger downward emission for structure B than that for structure A at wavelength 520 nm. However, as shown in Fig. 4.2, it is worth noting that the peak value of downward emission and the peak enhancement factor of downward emission for structure B at wavelength 520 nm are still smaller than those for structure A at wavelength 435 nm.

The reason is that the SPP resonant wavelength for structure B is still at wavelength 435 nm but its LSP resonant wavelength now occurs at 520 nm. Therefore, the resonant coupling effect between SPP and LSP at a common wavelength is impossible. Shown in Figs. 4.6 (a) and (b) are the distributions of the absolute values of magnetic field (a) and electric field (b) for structure D at wavelength 435 nm. The parameters for structure D are: r = 20 nm, d = 68 nm, t = 0.8 nm, and h = 20.8 nm. Note that the radius of each circular cylinder is purposely enlarged to see its effect. The peak value of downward emission for structure D is only 0.0007. This reduction results from the fact that the SPP effect is weak due to the large

(42)

distance from dipole to flat metal and the LSP effect is also weak because the LSP resonant wavelength now is shifted away from 435 nm.

In summary, we have synthesized a metallic nanostructure (structure A), for which the SPP and LSP can couple effectively, to enhance the downward emission of a nearby dipole. The peak enhancement factor of downward emission for three structures A, B and C is listed Table 4.1, from which several conclusions can be drawn as follows. The emission enhancement is best when both SPP and LSP are excited and coupled at a common resonant wavelength. It becomes worse when only one of SPP or LSP is effectively excited and hence the coupling between SPP and LSP is weak. The emission enhancement is worst when both SPP and LSP are weak.

4.2 Metallic nanostructures consisting of a half-space and separate circular cylinders

For further investigating the SP effect, we try to synthesize the metallic nanostructures consisting of a Ag half-space and two separate circular cylinders, as the insert of Fig. 4.7. In this section, we try to synthesize a metallic nanostructure for application to white-light LEDs.

Recently, most of the white-light LEDs in the market emit white light by

(43)

mixing blue light (wavelength 450 nm) and yellow light (wavelength 570 nm) [36]. Therefore, in this section the cost function is composed of downward emission at two wavelengths 450 nm and 570 nm.

In chapter 2 it is noted that the LSP resonant wavelength for a single thin Ag circular cylinder embedded in GaN is at 430 nm, near the desired blue light wavelength 450 nm. According to the assumed Drude model mentioned in chapter 2, the LSP resonant wavelength for a single thin Au circular cylinder embedded in GaN is at 540 nm, near the desired yellow light wavelength 570 nm. For the two separate circular cylinders considered here, the right one is Ag circular cylinder and the left is Au.

To evoke the SA process, we denote the downward emission by φ in all the following simulations, and define the cost function C in this section as C = -(α ×φ₄₅₀ + (1-α)×φ₅₇₀). (4.2) Here, φ₄₅₀ represents the downward emission at wavelength 450 nm, and φ570 represents that at wavelength 570 nm. For balancing the relative proportion between φ₄₅₀ and φ₅₇₀ in the cost function, we adjust parameter α to let α × φ₄₅₀ equal (1-α )× φ₅₇₀ during the iteration process. After an iteration, if the cost function of the new structure is less than the old one, we accept the new structure. Otherwise, we accept or

(44)

reject the new structure according to the SA criterion described in chapter 3. Therefore, if the new structure is accepted, both φ₄₅₀ and φ₅₇₀ are changed and parameter α has to be adjusted again. All the adjustable parameters and the way of changing the structure are described in section 3.2.2. Note that the vertical position of the dipole source is kept separate from the bottom of the circular cylinders by at least 5 nm.

Figure 4.7 shows the variation of cost function in the iteration process. Here the solid curve represents the cost function and the dashed curve represents the chosen temperature distribution used in the SA process. Due to the sensitivity and importance of initial state for convergence, firstly, we do a SA process with an initial state and choose a better structure from those ever appearing in the iterations as a new initial state. Next, we re-start the SA process with the new initial state for reducing the number of iterations. Figure 4.7 is the result from the secondary SA process. In the figure, we can find that the SA process indeed converges to the optimal structure. The final optimal structure labeled as structure E is obtained through the SA with r = 10 nm, d = 32 nm, t = 24 nm, and h = 40 nm. Note that the vertical distance from dipole to the bottom of circular cylinders is 6 nm, slightly larger than 5 nm.

(45)

Figure 4.8 shows the variations of downward emission φ₄₅₀ and φ₅₇₀ at two wavelengths in the iteration process. Here the yellow solid curve represents the downward emission for yellow light, and the blue dashed curve represents that for blue light. In the figure, we can find that both φ450 and φ₅₇₀ increase in the iteration process and reach the steady states.

Figure 4.9 shows the spectra of total emission and downward emission for structure E. The total emission is defined as the total outgoing power through a square of 2.5 nm in dimension centered at the dipole position in all the following simulations. Here the blue dashed curve represents the total emission, and the green solid curve represents the downward emission. In Fig. 4.9, we can find that there are two peaks in the spectrum of downward emission. One is at wavelength 450 nm with φ₄₅₀= 0.0072, and the other is at wavelength 565 nm with φ₅₆₅= 0.0181. This result is very consistent with our goal for SA.

Figure 4.10 shows the spectra of enhancement factor of total emission and downward emission for structure E. The enhancement factor of downward emission in this section is defined as the downward emission across the surface of y = -(60 + h) nm divided by the similar downward emission when the whole structure in the insert of Fig. 4.7 is