金屬孔洞陣列造成增強性兆赫輻射穿透時孔洞內介質的效應

(1)

國立交通大學

光電工程研究所

碩士論文

金屬孔洞陣列造成增強性兆赫輻射穿透時

孔洞內介質的效應

Effects of Hole material on Enhanced Terahertz

Transmission through Metallic Hole Arrays

研究生：羅誠

指導老師：潘犀靈教授

(2)

金屬孔洞陣列造成增強性兆赫輻射

穿透時孔洞內介質的效應

Effects of Hole material on Enhanced Terahertz

Transmission through Metallic Hole Arrays

研究生：羅誠 Student: Cheng Lo

指導老師：潘犀靈教授 Advisor: Prof. Ci-Ling Pan

國立交通大學光電工程研究所

碩士論文

A Thesis

Submitted to Department of Photonics & Institute of Electro-Optical Engineering

College of Electrical Engineering National Chiao Tung University In partial Fulfillment of the Requirements

for the Degree of Master of Engineering

In

Electro-Optical Engineering

July 2005

Hsinchu, Taiwan, Republic of China

(3)

(4)

摘要

二維週期次波長孔洞陣列的金屬薄板(2D-MHA)可以在某些特定的兆赫頻率(THz)造成異常高的穿透，這種增強性的穿透被認為是由於入射光與在金屬與介電質之間的表面電漿極化子(SPP)發生了共振耦合。我們利用實驗與模擬研究 2D-MHA 中洞裡介質在穿透特性中的角色，發現當洞裡填入介電質材料時有新的現象產生，當洞裡填入 UV 膠時，兆赫波的穿透特徵與 2D-MHA 的厚度有強烈的相關，且改變 2D-MHA 週遭介質的穿透特質也和以往不同；當孔洞深度較深時， SPP 展現出兩種不同的模態 (耦合與非耦合)，這也在利用時域有限差分法(FDTD)的電場模擬結果中獲得確認。

(5)

Abstract

Metal plates perforated with two-dimensional periodic arrays of sub-wavelength holes (2D-MHAs) exhibit high transmission at selective terahertz frequencies. This mechanism of the enhanced transmission characteristic is usually attributed to the resonant coupling between the incident light and surface plasmon-polaritons (SPPs) which are located at the metal-dielectric surfaces. We experimentally and numerically investigate the role of material in the holes on transmission characteristics of the 2D-MHA. New phenomena, however, appear when holes of the 2D-MHA are filled with dielectric material. We find distinctive THz transmission characteristics depending on the thickness and properties of the adjacent medium for the 2D-MHA of which holes are filled with UV-gel. For deeper holes, SPP can exhibit two distinct modes, coupled and uncoupled types as confirmed in the simulated results for the electric field using finite-difference time-domain (FDTD) algorithm.

(6)

Acknowledgement (誌謝)

碩士班兩年的日子真的很快，回顧這兩年，真的讓我成長了很多。懵懵懂懂地踏入了超快光電的領域，如今已經寫完了論文，今年似乎少了畢業的感傷，但對接下來的人生卻有著更多的期待。首先要感謝我的指導教授潘犀靈老師，從大學專題生時就一直受到他的照顧，口試委員賴暎杰老師、趙如蘋老師、洪勝富老師在口試時所給的建議也讓我獲益良多，光電所裡謝文峰老師、許根玉老師在問題上的教導也讓我的觀念更加清晰;實驗室裡的夥伴們:專題時期就帶我的藍玉屏學姐跟黃銘傑學長，學姊教導的光學儀器、雷射知識以及對光技巧讓我受用無窮，總是被我纏著一直問問題卻都耐心替我解答的老劉學長，你是我在 THz 遇到困難時最大的靠山，教我許多光學觀念的 chuck，提供良好雷射的 moya，教導我很多實驗室事務的 mika，帶領我進入 THz 領域的之揚，當初每次深夜時你的傾囊相授，讓我一步步地熟悉系統，是我碩士班最扎實的基礎，一起討論 THz 以及光子晶體卻常常找我集合的卓帆隊長，你是我心情不好時最好的傾訴對象，教導我怎麼製作液晶樣品的家任，提供我很多意見的昭遠、阿達、信穎、宇泰和已經畢業的上屆學長姊們，以及一起打拼衝刺，一起討論一起加油的同學們:跟我難兄難弟的宗翰、一起架設系統的阿隆、敎我很多電腦知識的小冷、很會消失但超聰明的小壯、憨厚正直的小高、查資料超強的 cc 和多才多藝的仔仔，還有可愛的學弟妹: 乃今、國騰與佳瑩，是你們讓我在實驗室有著無數歡樂的回憶;電子所楊玉麟學長在 FullWAVE 軟體上的教導，機械工廠莊師傅在樣品製作上的幫忙，都是我論文後的最大推手。家人們的陪伴更是我最大的力量，媽媽總是無時無刻地掛念著我，回家時陪著我聊天談心，讓我可以在學業上無憂無慮地衝刺，是我最安定的力量，大姐、二姐、姊夫你們都是我最親愛的家人，當然還有詩敏，每天跟妳講電話以及每次的見面都是我最快樂的時光，這麼多年來的關懷與包容是我最溫暖的寶貝。接下來的人生或許更加辛苦，但我會帶著你們滿滿的祝福開心地走下去﹗ 羅誠 2005 年 7 月于新竹交大

(7)

List of Figures

Fig. 1-1 Electromagnetic spectrum. The terahertz region ranges from frequencies of about 100 GHz to 10 THz. Refer to http://www.wse.jhu.edu/~cmsd/Thz/

Fig. 2-1 Simulated THz time domain waveform by solving the differential equations in Drude-Lorentz model referring to [13].

Fig. 2-2 Schematic diagram of a THz-TDS spectrometer using a femtosecond laser and photoconductive antennas to generate and detect THz waves.

Fig. 2-3 Typical current response J(t) of a photoconductive antenna to a short optical excitation pulse.

Fig. 2-4 Principle of photoconductive sampling. The photoconductive switch acts as a sampling gate that measures the waveform voltage V(t) within the sampling time . Fig. 2-5 Working principle of THz-TDS. By changing the optical delay between the optical pulse triggering the sampling gate and the waveform, the entire waveform can be mapped out sequentially in time.

Fig. 2-6 Schematic diagram of the 2D-MHA. There are three main parameters, such as hole diameter d, lattice constant s, and thickness of the metal plate t.

Fig. 2-7 A p-wave propagates along the interface between dielectric material ( )ε₁ and metal ( )ε₂ in the x-direction when z>0.

Fig. 2-8 The dispersion relation of the incident electromagnetic wave and the surface plasmon. Medium1: the dispersion curve of EM waves propagate in dielectric material; SP: the dispersion curve of surface plasmon.

Fig. 2-9 Yee’s mesh

Fig. 3-1 Terahertz Time-Domain Spectroscopy in our lab.

Fig. 3-2 The photo of our THz-TDS system. The photoconductive emitter and detector forms a symmetry type can be seen.

Fig. 3-3 Three different silicon design. (a) the non-focusing hemispherical design. (b) the hyper-hemispherical lens, and (c) , the collimating substrate lens.

Fig. 3-4 The schematic diagram of silicon lens and GaAs substrate. The radius of silicon lens is 6.75mm and the total thickness is 8.35mm. The thickness of the wafer is 0.355mm. The imaging point of the dipole antenna through the silicon lens is located in the left side with L distance.

Fig. 3-5 Photo of the setup of the THz emitter and the object lens. The emitter is combined with the object lens in one stage.

Fig. 3-6 Schematic diagram of the transmission properties of THz wave through MHA. THz wave is parallel between two parabolic mirrors and normal incident to the sample.

(10)

and thickness of the metal plate t=0.5 mm Fig. 3-8 The photo of the real MHA sample.

Fig. 3-9 (a) Procedure 1 of fabricating the MHA which holes filled with UV-gel. Let the MHA lay on the UV-gel from top to bottom, and then a MHA filled with UV-gel can be obtained.

Fig. 3-9 (b) Procedure 2 of fabricating the MHA which holes filled with UV-gel. Devising a mold for fixing the MHA.

Fig. 3-9 (c) Procedure 2 of fabricating the MHA which holes filled with UV-gel. Using the miller to drill the MHA, and to obtain a thickness t we want.

Fig. 3-10 The photo of the MHA with expected thickness which holes filled with UV-gel.

Fig. 3-11 (a) Dielectric function of Aluminum versus frequency. Fig. 3-11 (b) Refractive index of Aluminum versus frequency. Fig. 4-1 Free space THz time domain waveform

Fig. 4-2 THz frequency domain spectrum

Fig. 4-3 (a) The time domain signal compared with the reference Fig 4-3. (b) The frequency domain signals via the FFT from (a)

Fig. 4-3 (c) The power transmittance of the MHA. A obvious characteristic of band-pass filter, and the magnitude of the peak almost get up to 100%.

Fig. 4-4 (a) Simulated result: The time domain signal compared with the reference Fig. 4-4 (b) Simulated result: The frequency domain signals via the FFT from (a) Fig. 4-4 (c) Simulated result: The power transmittance of the MHA. The spectral linewidth of the peak is wider slightly than the experiment, but the peak frequency is also located at 0.3THz nearby.

Fig. 4-5 (a) Using SPP model to estimate the transmission peak in the triangular MHA. The peak frequency calculated from eq.(2.2.26) is 0.348THz right to the observed peak at 0.301THz

Fig. 4-5 (b) Using SPP model to estimate the transmission peak in the cubic MHA. The peak frequencies calculated both are right to the the observed peaks.

Fig. 4-6 (a) Altering the thickness of the first sample. The structure of the first sample is s =1.13mm, d =0.68mm.

Fig. 4-6 (b) Altering the thickness of the second sample. The structure of the second sample is s =0.99mm, d =0.56mm.

Fig. 4-7 (a) The dispersion relation of the UV-gel in THz region. The real part of refractive index is about 1.68 with almost non-dispersion.

Fig. 4-7 (b) The imaginary part of refractive index of the UV-gel in THz region. The attenuation seems to be negligible in our case.

(11)

Fig. 4-9 Power transmittance of the MHA which holes filled with UV-gel. The transmission peak broadened and multi-peak features are observed.

Fig. 4-10 Altering the thickness of MHA when the holes filled with UV-gel. All peaks shift to high frequencies and disappear when they approach the diffraction limit. Fig. 4-11 FDTD simulated results: Altering the thickness of MHA when the holes filled with UV-gel. It shows the same trend with experimental results in Fig. 4-9 Fig. 4-12 (a) The first peak frequency as the function of the thickness

Fig. 4-12 (b) The second peak frequency as the function of the thickness Fig. 4-12 (c) The first valley frequency as the function of the thickness

Fig. 4-12 (d) The spacing between the first and the second peaks as the function of the thickness

Fig. 4-13 (a) 100µm-thick MHA which holes filled with UV-gel attach different layers of ScotchTM tapes on the incident side. The number of tapes is from zero to five.

Fig. 4-13 (b) 100µm-thick MHA which holes filled with UV-gel attach different layers of ScotchTM tapes on the incident side. The number of tapes is four, six, eight and ten.

Fig. 4-13 (b) 100µm-thick MHA which holes filled with UV-gel attach different layers of ScotchTM tapes on the incident side. The number of tapes is from thirteen to fifteen.

Fig. 4-14 (a) 400µm-thick MHA which holes filled with UV-gel attach different layers of ScotchTM tapes on the incident side. The number of tapes is from zero to five.

Fig. 4-14 (b) 400µm-thick MHA which holes filled with UV-gel attach different layers of ScotchTM tapes on the incident side. The number of tapes is from six to ten. Fig. 4-15 (a) Simulated z component Ez of the electric field amplitude for 100µm-thick MHA filling with UV-gel. The incident CW wave at 0.191THz.

Fig. 4-15 (b) Simulated z component Ez of the electric field amplitude for 100µm-thick MHA filling with UV-gel. The incident CW wave at 0.276THz.

Fig. 4-16 (a) Simulated Ex and Ez for the 400µm-thick 2D-MHA at 0.191THz. (1)-(6) shows a cycle of propagation.

Fig. 4-16 (b) Simulated Ex and Ez for the 400µm-thick 2D-MHA at 0.276THz. (1)-(6) shows a cycle of propagation.

(12)

1 Introduction

Every region of frequency domain is a treasury in nature. The terahertz region of the electromagnetic spectrum ranges from frequencies of about 100 GHz to 10 THz (10 x 1012 Hz) as shown in Fig. 1-1 has wide and unexplored applications. This corresponds to wavelengths between about 3 and 0.03 mm, and lies between the microwave and infrared regions of the spectrum. At lower frequencies, microwaves and millimeter-waves can be generated by "electronic" devices such as those components in mobile phones. At higher frequencies, near-infrared or visible light is generated by "optical" devices such as semiconductor lasers, in which electrons emit light when they jump across the semiconductor band gap. Unfortunately, neither electronic nor optical devices cannot conveniently make work in the terahertz region because the terahertz frequency range sits between the electronic and optical regions of the electromagnetic spectrum.

Fig. 1-1 Electromagnetic spectrum. The terahertz region ranges from frequencies of about 100 GHz to 10 THz.

(13)

This situation changed in the 80's with the appearance of ultrashort pulse lasers of about 100 femtosecond pulse duration. Ultrashort pulses can generate broadband THz radiation when they impinge on semiconductors. Mourou and Auston et al. first demonstrated generation and detection of pulsed THz radiation by photoconducting switch with advantages of time resolution of picosecond and sensitivity enhanced by phase-lock technique [1-2]. In 1996, Zhang et al. developed free-space electro-optic sampling (FS-EOS) technique to enhance signal to noise ratio (S/N ratio) up to 10000 and to achieve much large dynamic range [3]. These techniques have now developed to a level for spectroscopy and sensing. The spectroscopic technique using pulsed THz radiation is called “THz time-domain spectroscopy (THz-TDS)” [4].

On the other hand, microstructured devices, such as dichroic filters [5], terahertz plasmonic filters [6], and terahertz photonic crystals [7] have received growing interests. These frequency-selective components play an important role in the development of terahertz technology. In 1998, Ebbesen et al. first explored extraordinary transmission of light in optical frequencies through a metal film perforated by an array of sub-wavelength holes (MHA) [8]. Many researchers have paid much great concentration on this field. Since the period of hole arrays must lie on the order of wavelength, so the most difficult challenge to studying periodical hole array is the manufacture in visible light. Fortunately, the wavelength in THz region is much longer than visible light. For this reason, it is much easier to fabricate periodical hole arrays in THz region. Furthermore, the loss of metal in THz region is negligible so metals are

(14)

However, theories to describe this phenomenon are still ambiguous. The most scientists consider these unusual phenomena are caused by the coupling between the incident waves and the surface plasmon polaritons (SPPs) [10]. Base on the theory of SPPs, the maximum transmittance peak should to be independent of the hole’s diameter and the material stuffed with the holes [11]. However, in our experimental results, when the holes filled with the UV-gel, the number of the transmittance peak is changed from single to multeity, and the peak frequency will change with the thickness of the samples. Further, we attach translucent ScotchTM tapes on the incident side of MHA filled with UV-gel. Peaks shift to the left and decrease as the SPP resonances approach the cutoff frequencies. These phenomena are analogous with previous reports elucidated in terms of the SPP [12]. However, upon increasing the number of the tape, a side peak on the high frequency side grew gradually and red-shifted. For the number of tapes up to fifteen layers, the trend of shift still persists. We also discovered that different peaks have different SPP-like surface wave modes as confirmed in the simulated results using Finite-Difference Time-Domain (FDTD) for the electric field. These findings are important for understanding the fundamental mechanisms of enhanced transmission and will be vital in designing devices based on this effect.

(15)

2 Basic Theories and Conceptions

2.1 Terahertz (THz) Field

Terahertz (THz) fields is a generic term for electromagnetic waves

with a spectrum between 0.1 and 10 THz (where 1 THz is 1012

cycles/second). From the viewpoint of wavelength, it includes millimeter and submillimeter waves. THz signals were until recently an almost unexplored area of research due to the difficulties in generation and detection of electromagnetic fields at these wavelengths. The use of ultrafast lasers to generate subpicosecond pulses of electromagnetic radiation, THz pulses, has evolved into a very active research field during the past two decades. In the other hand, THz detection is also a quite difficult work since the power of emitted THz signals is weak. Owning to the high signal-to-noise ratio, the use of photoconductive antennas to detect THz wave is a very popular method.

2.1.1 Generation of THz Radiation Using Photoconductive Antennas

When the photoconductive antenna is illuminated by ultrashort optical pulse, where the photon energy is greater than the bandgap of the semiconductor, a planar photoconductor absorbs the incident light, coherently exciting free electron-hole pairs. Photo carrier acceleration results from an applied bias on the electrodes, which produces a transient surface current. The expression can be described from the current-surge model [13].

(16)

Current-surge model:

Initially, the radiating source is defined as time-varying parameter, including charge densityρ

(

x,y,z,t

)

, current density JK

(

x,y,z,t

)

, electric

fieldEG

(

x,y,z,t

)

, or magnetic flux BK

(

x,y,z,t

)

. And then, it is necessary to

construct Maxwell’s equation before deducing current-surge model. Maxwell’s equation: t B E ∂ ∂ − = × ∇ K K (Faraday’s Law) (2.1.1) E ρ ε ∇ ⋅ =K (Gauss’ Law) (2.1.2) t D J H ∂ ∂ + = × ∇ K K K (Ampere’s Law) (2.1.3) 0 B ∇ ⋅ =K (2.1.4) From (2.1.1) and A BK =∇× K (2.1.5)

(

)

_⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ − × ∇ = × ∇ ∂ ∂ − = ∂ ∂ − = × ∇ t A A t t B E K K K K 0 = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + × ∇ ⇒ t A E K K (2.1.6) And then, some non-vector value, V, is induced in the equation.

From (2.1.6) Set t A E V ∂ ∂ + = ∇ − K K (2.1.7) And then,∇×

(

−∇AK

)

=0 (2.1.8) From (2.1.7), t A V E ∂ ∂ − −∇ = K K (2.1.9) Next, the two inhomogeneous wave equations written in terms of A

(17)

and V could also be deduced from the inhomogeneous Maxwell equations: From (2.1.3), E D B H t E J B t D J H K K K K K K K K K K ε µ ε µ = = ∂ ∂ + = × ∇ ⇒ ∂ ∂ + = × ∇ and As (2.1.10)

From (2.1.10), (2.1.5) and (2.1.9) can be written as

(

)

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ − ∇ − ∂ ∂ + = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + = × ∇ × ∇ t A V t J t E J A K K K K K ε µ ε µ (2.1.11)

(

)

2 2 2 2 2 2 V A A A J t t A V A J A t t µ µε µε µε µ µε ∂ ∂ ⎛ ⎞ ∇ ∇ ⋅ − ∇ = − ∇_⎜ _⎟− ∂ ∂ ⎝ ⎠ ∂ ⎛ ∂ ⎞ ⇒ ∇ − = − + ∇ ∇ ⋅ +_⎜ _⎟ ∂ ⎝ ∂ ⎠ K K K K K K K K (2.1.12) From (2.1.2),

( )

(

)

2 D A E V t V A t ρ ε ε ρ ρ ε ∇ ⋅ = ⎡ ⎛ ∂ ⎞⎤ ⇒ ∇ ⋅ = −∇ ⋅_⎢ _⎜∇ + _⎟_⎥= ∂ ⎝ ⎠ ⎣ ⎦ ∂ ⇒ ∇ + ∇ ⋅ = − ∂ K K K K (2.1.13) Set A V 0 t εµ∂ ∇ ⋅ + = ∂ K (Lorentz gauge) (2.1.14) So that, (2.1.12) becomes as J t A A K K K µ µε =− ∂ ∂ − ∇2 2₂ _(2.1.15) From (2.1.14), A V t µε ∂ ∇ ⋅ = − ∂ K

As a result, (2.1.13) can be written as

ε ρ µε =− ∂ ∂ − ∇2 2₂ t V V (2.1.16)

(18)

and (2.1.16) are demonstrated. These equations here are used to determine a functional, time dependent form of the radiated electric field in the far field.

From (2.1.3), the continuity equation of the free carriers, which is generated in the biased semiconductor after the absorption of an optical pulse, is obtained.

(

H

)

J D J 0 t t ρ ⎛ ∂ ⎞ ∂ ∇ ⋅ ∇× = ∇ ⋅_⎜ + _⎟= ∇ ⋅ + = ∂ ∂ ⎝ ⎠ K K K K (2.1.17) Actually, the current in the bias photoconductor is a transverse current, which is parallel to the surface of the photoconductor and perpendicular to the direction of propagation, so that

0

J

∇ ⋅ =K (2.1.18) Equation (2.1.17) and (2.1.9) imply that the charge density dose not vary in time and not contribute to the time-dependent radiated electric field. As a result, from (2.1.9) the electric field is

( )

A

( )

t t t EK_rad K ∂ ∂ − = (2.1.19) The solution to the wave equation (2.1.15) and hence for the vector potential AK leads to the expression for the time-dependent radiated

electric field EK_rad =( trK, ) at a displacement rK from the center of the

photoconductor:

( )

da r r c r r t r J t c t r E s rad ₋ _′ ′ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − ′ − ′ ∂ ∂ − = _∫ _K _K K K K K K , 4 1 , ₂ 0 πε (2.1.20) where ε₀ is the permittivity of free space, c is the speed of light in

(19)

vacuum, JK_s is the surface current in the photoconductor evaluated at the

retarded time, and d ′a is the increment of surface area at a displacement

rK′from the center of the emitter. Integration is taken over the optically

illuminated area of the photoconductor. In the far field,

r r r n r r r ⎟≈ ⎠ ⎞ ⎜ ⎝ ⎛ ₋ ⋅ ′ = ′ −K K K ˆ 1 (2.1.21)

At the same time, the gap between the electrodes of the photoconductor is assumed to be uniformly illuminated by the optical Therefore, the surface current JK_s can be assumed to be constant at all

points on the surface of the emitter. And then, the radiated electric field can be written as

( )

(

)

⎟⎠ ⎞ ⎜ ⎝ ⎛ − + + − = c r t J dt d z y x A c t r EK_rad K ₁_/₂ K_s 2 2 2 2 0 4 1 , πε , (2.1.22)

where A is the illuminated area of the emitter. It is considered that the radiation emitted (and detect) on axis (i.e. x=y=0), and let ⎟

⎠ ⎞ ⎜ ⎝ ⎛ − → c z t t . Thus,

( )

( ) 4 1 , ₂ 0 t J dt d z A c t r EK_rad K K_s πε − ≅ (2.1.23) There is a result for the surface current

( )

_{( )}

1 1 0 ₊ + = n t E t J s b s s σ η σ K K (2.1.24) where n is the index of refraction of the photoconductor at submillimeters,

0

η is the impedance of free space which is equal to 377 , and the surface conductivity defined as

Ω

(

)

∫ ′ ′ ′ ⎢⎡− − ′ ⎥⎤ − =e R t dtm t t I t t t τ ω σ (1 ) ( , ) ( )exp (2.1.25)

(20)

where e is the electric charge, R is the optical reflectivity of the photoconductor, =ω is the photon energy , m( tt, ′) is the carrier

mobility at time t of a carrier created at time t′ , is the

time-dependent optical intensity, and

opt

I

car

τ is the lifetime of the excited carriers. For the present derivation, a constant carrier mobility is assumed:

m t t

m( , ′)= (2.1.26)

Also, the carrier lifetime is long, i.e. τ_car →∞. Finally a Gaussian intensity profile of the optical beam is assumed:

⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛− ′ = ′) ₀exp ₂2 ( τ t I t I_opt . (2.1.27)

With these assumption s the surface conductivity becomes

(

)

∫ ∞ − ⎟⎟⎠ ⎞ ⎜⎜ ⎝ ⎛− ′ ′ − = t s t m t d I R e 2 2 0 exp 1 τ ω σ = (2.1.28)

Equations (2.1.23), (2.1.24), and (2.1.28) lead to the expression

( ) _/ 2 2 0 0 2 2 0 2 0 ) exp( ) 1 ( ) 1 ( 1 exp 1 4 ) ( − ∞ − ⎥⎦ ⎤ ⎢ ⎣ ⎡ − + − + × ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − − − = _∫ττ ω η τ ω πε t b rad x dx n m I R e t m I R e z A c E t E = = K K _(2.1.29)

As the comparison with the result from experiment, it is necessary to rewrite the equation in terms of the experimental parameters , the bias electric field applied across the photoconductor, and , the incident optical fluence, which is defined as

b E opt F

∫

−∞ ∞ ≡ = − = A E I dt t I F_opt

π

τ

opt

τ

2 0 2 0exp( ) (2.1.30)

(21)

where is the average optical energy and A is the area of the optical beam, where r is the measured

opt

E

e

1_{radius of the optical beam. Set}

2 2 0 (1 ) (m / ) 4 Ae R m B s J c z πε ω π − = = (2.1.31) ). / (m ) 1 ( ) 1 ( 2 0 J n m R e D π ω η = + − = (2.1.32) And then, the electric field in the far field can be written as

. ) exp( 1 ) exp( ) ( 2 2 2 2 − ∞ − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + × − − = _∫τ τ τ t opt opt b rad DF x dx t F BE t E (2.1.33)

We can also discuss on the viewpoint of the carrier dynamics in semiconductor to analyze the THz generation by Drude-Lorentz model when using the photoconductive antenna [14, 15].

Drude-Lorentz model:

When a biased semiconductor is pumped by a laser pulse with photon energies greater than the band gap of the semiconductor, electrons and holes are created in the conduction band and valence band, respectively. The time dependence of carrier density is given by the following equation:

c

dn n

G

dt = −

τ

+ (2.1.34) where n is the density of the carrier, G is the generation rate of the carrier by the laser pulse, and τc is the carrier trapping time. The

(22)

acceleration of electrons (holes) in the electric field is given by , , , , e h e h e h s e h dv v q E dt = − τ +m (2.1.35)

where is the average velocity of the carrier, is the charge of an electron (a hole), is the effective mass of the electron (hole),

, e h v q_{e h}_, , e h m s

τ is the momentum relaxation time, and E is the local electric field. The subscripts e and h represent electron and hole, respectively. The local electric field E is smaller than the applied bias electric field, , due to the screening effects of the space charges,

b E b

P

E

αε

=

−

(2.1.36)

where P is the polarization induced by the spatial separation of the electron and hole, ε is the dielectric constant of the substrate, and α is the geometrical factor of the photoconductive material. The geometrical factor α is equal to three for an isotropic dielectric material. It is noted that both, the free and trapped carriers contribute to the screening of the electric field. The time dependence of polarization P can be written as

r dP P

J

dt = −τ + (2.1.37)

where τr is the recombination time between an electron and hole.

In eq. (2.1.37), J is the density of the current contributed by an electron and hole,

h

J

=

env

−

env

_e (2.1.38) where e is the charge of a proton. The change of electric currents

(23)

leads to electromagnetic radiation according to Maxwell’s equations. In a simple Hertzian dipole theory, the far-field of the radiation is given by THz E THz

J

E

t

∂

∝

∂

(2.1.39)

To simplify the following calculations, we introduce a relative speed v between an electron and hole,

h

v

= −

v

_e (2.1.40) Then the electric field of radiation can be expressed as

THz

n

v

E

ev

en

t

∂

∝

+

∂

(2.1.41)

The first term on the right-hand side of eq. (2.1.41) represents the electromagnetic radiation due to the carrier density change, and the second term represents the electromagnetic radiation which is proportional to the acceleration of the carrier under the electric field.

The simulated THz waveform in time domain by solving the differential equations in Drude-Lorentz model is shown in Fig. 2-1 referring to [13]. The momentum relaxation time is 30 fs, the carrier trapping time for the emitter is 0.5 ps, and the carrier recombination time

r

τ is taken as 10 ps. The carrier generation density is 1018 cm-3 and the applied bias electric field is 20 KV/cm.

(24)

Fig. 2-1 Simulated THz time domain waveform by solving the differential equations in Drude-Lorentz model referring to [13].

2.1.2 Detection of THz Radiation Using Photoconductive Antennas

Actually, the measurement of THz signal is the measurement of the photocurrent induced in the detector by the simultaneous arrival of the THz pulse and the gating optical pulse [13]. Because of the non-instantaneous conductivity of the detector σd, a measurement of the

photocurrent j induced by the arrival of the peak of the radiated electric field of the THz pulse at the detector at time is actually a convolution of the radiated electric field and the conductivity of the detector after time

:

p

t

p

(25)

∫

∞ ′ − ′ ′ = p t p d rad p E t t t dt t j( ) ( )σ ( ) . (2.1.42)

For the detector used in the experiments, the conductivity can be assumed to have the form

⎪ ⎩ ⎪ ⎨ ⎧ ′ ≤ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡− ′− ′ > = − ′ t t t t t t d p d p d p 0 p t if ) ( exp t if 0 ) ( τ τ σ σ (2.1.43)

where σ₀ is the peak conductivity of the detector. The rise in the conductivity can be assumed to be rapid because of the gating optical pulse has a short duration. Therefore, it is assumed that τd >>0, and

Eq.( 2.1.43) can be written as

⎪⎩ ⎪ ⎨ ⎧ ′ ≤ ′ > = − ′ t t t t d p d p 0 p t if t if 0 ) ( τ σ σ (2.1.44) Substitution of Eq. (2.1.44) into Eq. (2.1.42) gives

∫

∞ ′ ′ = p t rad d p p t E t dt j ( ) 0 ( ) τ σ _(2.1.45)

If it is assumed that the radiated waveform of the THz pulse is symmetric, then the measured photocurrent becomes

∫

∞ ∞ − ′ ′ ≈ E t dt j _rad d ) ( 2 0 τ σ (2.1.46) Substitution of Erad(t′) as given by Eq.(2.1.41 )into Eq. (2.1.46) gives

t d dx x DF t F BE j t opt opt b d ′ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − + × ′ − − = − ∞ ∞ − ′ ∞ −

∫

2 / 2 2 2 0 _exp( ₎ ₁ _exp( ₎ 2 τ τ τ τ σ _(2.1.47)

However, because of the internal reflection the back side of the photoconductor, the amplitude of the radiated electric field is attenuated.

(26)

The Fresnel formulas for reflection at a plane interface and conservation of energy, which indicate that the amplitude of the transmitted pulse is decreased relative to its initial value by the amount

1 2

+

n

n _(2.1.48)

Evaluation of Eq. (2.1.47) and multiplication of the result by Eq. (2.1.48) lead to the following expression for the measured detector photocurrent:

0 / 1 ( / ) s b s F F j CE F F σ = − + (2.1.49) where , ) 1 ( ) 1 ( 0e R m n F_s − + = η ω = _(2.1.50) . 4 2 z c n A C τ η πε = 0 0 d (2.1.51) As shown above, n is the index of refraction of the emitter at submillimeter wavelengths, =ω is the photon energy, η₀ is the

impedance of free space, R is the optical reflectivity of the emitter, A is the illuminated area of the emitter, ε₀ is the permittivity of free space, c is the speed of light, z is the emitter-detector separation, and -τd is the

lifetime.

Eq. (2.1.49) can be rewritten as

0 / ( ) 1 ( / ) s rad b s F F j E peak CE F F σ = = − + (2.1.52)

where is the experimentally measured value associated with the peak of the radiated electric field and the right

side, , is the form of the incident excitation

optical fluence and the bias field applied to the emitter. The constant

) ( peak E_rad

(

)

[

{

/ s /1 ( / s) b F F F F CE + −

]

}

b E

(27)

C provides a measure of the fraction of the bias field extracted in the far field, whereas is the optical fluence required for extracting half the maximum field.

s

F

2.1.3 Terahertz Time-Domain Spectroscopy (THz-TDS)

Figure. 2-2 shows a schematic diagram of a THz-TDS spectrometer [4]. It consists of a femtosecond laser source (1), a beam splitter that divides the laser into two branches, a THz emitter (2), focusing and collimating optics (3), the sample (4), THz detector (5), a variable delay line (6) that alters the optical delay between the THz emitter and detector, and a lock-in amplifier (7). A computer (8) controls the delay line and displays the detector photocurrent versus path length.

The sources and detectors of THz-TDS consist of the same building blocks. Both are based on the photoconductive antenna, which consists of a semiconductor bridging the gap in a transmission line structure deposited on the semiconductor substrate. The response of the voltage-biased photoconductive antenna to a short optical pulse focused into the gap between the two metal contacts is illustrated in Fig. 2-3.

(28)

(1) Beam Splitter (2) (3) (5) (4) (6) (7) (8) (1)Femtosecond Laser (2)THz Emitter (3)Collimating Optics (4)Sample (5)THz Detector (6)Optical Delay Line (7)Lock-In Amplifier (8)Computer (1) Beam Splitter (2) (3) (5) (4) (6) (7) (8) (1)Femtosecond Laser (2)THz Emitter (3)Collimating Optics (4)Sample (5)THz Detector (6)Optical Delay Line (7)Lock-In Amplifier (8)Computer

Fig. 2-2 Schematic diagram of a THz-TDS spectrometer using a femtosecond laser and photoconductive antennas to generate and detect THz waves.

Optical

Pulse

Photocurrent

Time

( )

J t

Optical

Pulse

Photocurrent

Optical

Pulse

Time

( )

J t

Photocurrent

Fig. 2-3 Typical current response

J t

( )

of a photoconductive antenna to a short optical excitation pulse

.

(29)

The current through the switch rises very rapidly when injecting photocarries by the optical pulse, and then decay with a time constant related to the carrier life time of the semiconductor. The transient photocurrent radiates into free space according to Maxwell’s equations,

( )

J t

( )

J t

E t

t

∂

∝

_{∂ (2.1.53)}

Because of the time derivative, the radiated field is dominated by the rising edge of the photocurrent transient, which is much faster than the decay. Long tail of the photocurrent decay, which occur in most semiconductor without high defect density, are largely irrelevant to the radiated field.

To convert photoconductive antenna to a detector of short electrical pulses, a current-to-voltage amplifier is connected across the photoconductor, replacing the voltage bias. The electric field of the incident THz wave now provides the driving field for the photocarriers. Current flows through the switch only when both the THz wave and photocarriers are present. Because electronics are not fast enough to measure the THz transients directly, repetitive photoconductive sampling is used (Fig. 2-4). If the photocarrier lifetime,

τ

is much shorter than the THz pulse, the photoconductive switch acts as a sampling gate which samples the THz field within a time

τ

.

(30)

Time

τ

THz Waveform Sampling Gate

( )

V t

Time

τττ

THz Waveform Sampling Gate

( )

V t

Fig. 2-4 Principle of photoconductive sampling. The photoconductive switch acts as a sampling gate that measures the waveform voltage V(t) within the sampling time

τ

.

Because the laser pulses which launch the emitter and gate the detector originate from the same source, the photoconductive gate can be moved across the THz waveform by changing the optical delay line shown in Fig. 2-2 (6). Fig. 2-5 shows using femtosecond pulse to analyze the THz waveform whose duration is about picoseconds. We can get entire THz transient waveform by employing this technique without the need for fast electronic components.

(31)

Sampling Gate THz Waveform Delay Time Electric Field Sampling Gate Delay Time Electric Field THz Waveform

Fig. 2-5 Working principle of THz-TDS. By changing the optical delay between the optical pulse triggering the sampling gate and the waveform, the entire waveform can be mapped out sequentially in time.

2.1.4 Extraction of optical constant in THz-TDS

If the THz pulses before and after the sample insertion are denoted

as and , respectively, then the ratio of their Fourier

transforms can be given by the complex refractive index of the sample

ref E ( )t Esam( )t ( ) n ω as [16] 2 0

( )

( ( ) 1)

( ) ( ) exp(

)

( )

2 ( )

(( ( )) exp(

))

( ) exp( ( ))

sam as sa ref m l sa l

E

d n

t

i

E

c

dn

r

i

c

T

i

ω

_ω

ω

ω ω

ω

φ ω

=

−

=

×

=

∑

(2.1.54)

(32)

Since the field is a real quantity, the Fourier transform will be a complex Hermitian spectrum with an angular frequency

ref E ( )t ω , ( ) _real( ) _imag( ) E ω =E ω −iE ω , where E( )ω =E(−ω) . The t_as( )ω , t_sa( )ω , ( ) sa

r ω are the real Fresnel coefficients for amplitude transmission and

reflection at the sample surfaces. The d means the thickness of the sample. The exponential factor represents the phase shift due to propagation in the material of the sample, and the last term is the contribution of the multiple reflection. The integer m represents the multiple reflection concluded in the observational time domain. T( )ω and φ ω( ) are

experimentally obtained power transmittance and relative phase, respectively. Since it is difficult to obtain the complex refractive index

( ) ( ) ( )

n ω =n ω −iκ ω directly from (2.1.55), the expressions for the real and

imaginary parts of the complex refractive index of the sample are denoted as 1 0 2

( )

( ( )

arg( ( ) ( )

2 ( )

(( ( )) exp(

)) ))

m as sa l l sa

c

d

n

t

d

c

n

d

r

i

c

ω

φ ω

ω

ω ω

ω

+ =

=

+

−

t

∑

(2.1.56)

2 1 2 0

( )

ln(

)

2 _{2 ( )}

( ) ( )

m

(( ( )) exp(

))

l as sa sa l

c

T

d

_n

_d

t

r

i

c

ω

κ ω

ω

_{ω ω}

ω

+

ω

=

= −

−

∑

(33)

We use the value of the complex refractive index roughly estimated from the THz pulse in the time domain before and after the sample insertion as a starting point for an iterative loop wherein n( )ω and

are calculated in a self-consistent manner. By performing this cyclic procedure for only a few times, the value would lead to a convergence.

( )

(34)

2.2 Metallic Hole Arrays (MHAs)

Metal films with two-dimensional periodic arrays of subwavelength holes can exhibit extraordinary optical transmission characteristics [8]. This discovery attracts much attention because of its potential applications for subwavelength optics. The enhanced transmission is explained to be due to the resonant coupling of incident light with surface plasmon polaritons (SPPs) [10]. The enhanced peak was observed that it will be changed strongly with different geometry of the structure. There are three main parameters, such as hole diameter d, lattice constant s, and thickness of the metal plate t. The hole arrays can be cubic or triangular. Schematic diagram of two-dimensional metallic hole array (2D-MHA) with triangular lattice is shown below

s

t

d

(35)

Fig. 2-6 Schematic diagram of the 2D-MHA. There are three main parameters, such as hole diameter d, lattice constant s, and thickness of the metal plate t

Since every hole can be thought as a waveguide, there must exist a cut-off frequency in this sample. Moreover, when the frequency of the incident wave is larger than diffraction limit, it will be diffracted into the first lobe due to the periodic structure [5]. It will be discussed explicitly below.

2.2.1 Cutoff Frequency and Diffraction Limit

Consider a circular waveguide with metallic surface and radius a. Equations of time-harmonic electric and magnetic field can be written as following [17]: 2 2 2 2

0

0 E

k E

H

k H

∇

+

=

∇

+

=

(2.2.1)

For a straight and uniform circular waveguide, it is convenient to decompose the 3-D Laplacian operator ▽2 into two parts: ▽rf2 and ▽z2

for transverse and longitudinal components, respectively. For TM waves, Hz = 0 and Ez ≠ 0, all fields can be expressed in terms of Ez E ez0 z

γ −

= ,

where 0

z

E satisfies the homogeneous Helmholtz’s equation 2 0 2 2 0 2 0 2 0 ( ) 0 r z z r z z E k E or E h E φ φ

γ

∇ + + = ∇ + = 0 . (2.2.2) For TE waves, Ez = 0 and Hz ≠ 0, all fields can be expressed in terms of

0 z z z

(36)

equation listed above.

For TM waves in circular waveguides,

0

( , , ) ( , ) z

z z

E r

φ

z =E r

φ

e−γ . (2.2.3)

The solutions are

0 0 0 2 0 2 0 0 ( )cos '( )cos ( )sin ( )sin '( )cos 0 z n n r n n n n r n n n n z E C J hr n j E C J hr n h j n E C J hr n h r j n H C J hr h r j H C J hr h H φ φ φ β n n φ β _φ ωε _φ ωε _φ = = − = − = − = − = (2.2.4)

where jb is equal to g, Jn is a Bessel function, and Cn is a coefficient. The

eigenvalues of TM modes are determined from the boundary condition that 0 z E must vanish at r = a: ( ) 0 n J ha = (2.2.5)

For the lowest TM mode, TM01, 01

2.405 ( )h _TM

a

= , (2.2.6) which yields the lowest cutoff frequency:

01 01 ( ) _0.383 ( ) 2 TM c TM h f a

π µε

µε

= = (Hz) (2.2.7) For TE waves in circular waveguides,

0

( , , ) ( , ) z

z z

H r

φ

z =H r

φ

e−γ (2.2.8)

(37)

0 0 0 2 0 2 0 0 ' ( )cos ' '( )cos ' ( )sin ' ( )sin ' '( )cos 0 z n n r n n n n r n n n n z H C J hr n j H C J hr n h j n H C J hr n h r j n E C J hr n h r j E C J hr n h E φ φ

φ

β

_φ

β

_φ

ωµ

_φ

ωµ

_φ

= = − = = = = (2.2.9)

where Cn’ is a coefficient. The eigenvalues of TE modes are determined

from the boundary condition that the normal derivative of 0

z H must vanish at r = a: '( ) 0 n J ha = (2.2.10)

For the lowest TE mode, TE11, 11

1.841 ( )h _TE

a

= , (2.2.11) which yields the lowest cutoff frequency:

01 11 ( ) _0.293 ( ) 2 TE c TE h f a

π µε

µε

= = (Hz) (2.2.12) In conclusion, the cutoff frequency for an infinitely long cylindrical waveguide is defined by the hole diameter d

1.841 cutoff c d

ν

π

= _(2.2.12) where is speed of light in free space. c

In addition to the cutoff frequency, owing to the triangular hole array represents a 2-D grating. When electromagnetic waves with frequency

(38)

above 2 3 diff c s ν = (2.2.13) it will be diffracted into the first lobe. where is the hole spacing. s

2.2.2 Theories of Surface Plasmons (SPs)

In order to understand the characteristics of surface plasmons, we consider a p-wave propagates along the interface between dielectric material ( )ε₁ and metal ( )ε₂ in the x-direction. The metal locate in the area where z<0, such as in Fig. 2-7:

z

x

k

_x1

k

_x2

k

_z1

k

_z2 (Dielectric Material) (Metal) 1

( )

ε

2

( )

ε

Fig. 2-7 A p-wave propagates along the interface between dielectric material ( )ε₁ and metal ( )ε₂ in the x-direction when z>0.

(39)

1 1 1 1 1 1 1 1 1

(0,

,0) exp (

)

(

,0,

) exp (

y x z x z x z

)

H

i k x

k z

t

E

i k x

k z

t

ω

=

+

=

+

−

)

(2.2.14) when z<0 2 2 1 1 2 2 2 2 2

(0,

,0) exp (

)

(

,0,

) exp (

y x z x z x z

H

i k x k z

t

E

i k x k z

ω

t

ω

=

−

=

−

(2.2.15)

The above equations must satisfy the Maxwell equations,

1 1 0 0 1, 2 i i i i i i i i H E E H E H c t c t and i µ ∂ ε ∂ ∇ ⋅ = ∇ ⋅ = ∇ × = − ∇ × = ∂ ∂ = K K K K K K K K K K (2.2.16)

Then, using the boundary conditions,

1 2 1 2 1 1 2 2

x x y y z z

E

=

E

H

=

H

ε

E

=

ε

E

(2.2.17) We can obtain the following results:

2 ,

1 ,

]

)

(

[

)

(

2 1 2 2 2 1 2 1 2 1 2 1

=

−

=

+

=

≡

=

i

k

c

k

c

k

x i zi x x x

ω

ε

ω

(2.2.18)

This is the relative formula for the surface plasmons in metal. We can find that 2 1 1 2 1 2 1 2 1 ₎ (

ω

ε

ω

c c k_x > + = _(2.2.19) That is to say the wavevector of electric field in the surface plasmon is larger than the incident electromagnetic wave. Fig. 2-8 shows the dispersion relation of the incident electromagnetic wave and the surface plasmon.

(40)

x

k

c

=

1

ε

ω

l : ω k_x SP : x k + c = 2 1 2 1 ε ε ε ε ω 1 ε 1 ω + p x

k

c

=

1

ε

ω

l : ω k_x SP : x k + c = 2 1 2 1 ε ε ε ε ω 1 ε 1 ω + p

ω

x

k

1 1: c _x medium

ω

k

ε

= 1 2 1 2 : _x SP

ω

c

ε ε

k

ε ε

+ = x

k

c

=

1

ε

ω

k

c

=

1

ε

ω

k

c

=

1

ε

ω

k

c

=

1

ε

ω

l : ω k_x SP : x k + c = 2 1 2 1 ε ε ε ε ω 1 ε 1 ω + p

ω

x

k

1 1: c _x medium

ω

k

ε

= 1 2 1 2 : _x SP

ω

c

ε ε

k

ε ε

+ =

Fig. 2-8 The dispersion relation of the incident electromagnetic wave and the surface plasmon. Medium1: the dispersion curve of EM waves propagate in dielectric material; SP: the dispersion curve of surface plasmon.

Since the dispersion curve of surface plasmon always lies to the right of the incident electromagnetic wave so the usual incident photons can’t excite surface plasmon polaritons (SPPs). However, periodical holes in the metallic film can lead to the in-plane momentum which aids the photon energy coupled to surface plasmon polariton modes [18]. The corresponding relation for conservation of momentum is

sp x x

k

K

=

k

K

+

iG

K

+

jG

K

_y _(2.2.20)

where kK_sp is the wave vector associated with the SPP, is the transverse wave vector component of the incident radiation, and G

x

kK

x

GK K_y

are the wave vector components associated with the two-dimensional array, and i and j are integers. k_x =(2 / sin )π λ θ is the component of the

(41)

wave vector of the incident light that lies in the plane ( θ =0 for normal incidence), From the above calculation,

1 2

(

d m sp d m

k

c

ε ε

ω

ε

=

+

)

(2.2.21) where ω is the frequency of the incident electromagnetic radiation, is the speed of light in vacuum,

c

d

ε is the dielectric constant of the dielectric interface medium and εm is the (complex) dielectric constant

of the metal. This last quantity may be expressed as εm=εmr+iεmi ,

where εmr and εmi are the real and imaginary components of the

dielectric constant of the metal, respectively. The complex propagation constant for the SPP wave can be expressed as ksp =kspr +ikspi, where the

two individual components can be written as [19]

1 1 2 4 2 _{2 2} 1 2 2 2 2 1 2 1 1 2 2 2 4 2 _{2 2} ₂

(

)

[

] [

]

(

)

2 [

]

(

)

{2[

(

) ]}

d e e d mi spr sp mr d mi d d mi spi mr d mi e e d mi

k

n

c

k

c

ε

ε ε

ω

ε ε

ε

ε ε

ω

ε ε

ε

_ε

_{ε ε}

+ +

+

=

+

=

+

(2.2.22) with 2 2 2 e mr mi d mr

ε =ε +ε +ε ε . In the limit where εmi > εmr and

,

mr mi d

ε ε ε , which is true for most metals at THz frequencies, the real and imaginary components of the SPP wave vector can be greatly simplified. The real component of the SPP propagation constant may be approximated as

(42)

spr sp d

k

n

c

ω

_{ω ε}

=

_(2.2.23)

and the imaginary component of the SPP propagation constant may be approximated as 3 2

2

d spi mi

k

c

ε

ω

ε

=

_(2.2.24)

At normal incidence, Eq. (2.2.23) gives the position of the maxima:

1

2 2 ₂

max

a i

0

(

j

)

d

λ

=

+

−

ε

_(2.2.25) for a square lattice. For a triangular lattice a similar derivation yields

1 2 2 ₂ max 0

4 [ (

)]

3

d

a

i

ij

j

λ

=

+ +

−

ε

_(2.2.26)

(43)

2.3 Finite-Difference Time-Domain(FDTD) Algorithm

Maxwell’s partial differential equations of electrodynamics, formulated circa 1870, represent a fundamental unification of electric and magnetic fields predicting electromagnetic wave phenomena which Nobel Laureate Richard Feynman has called the most outstanding achievement of 19th-centry science. Now engineers and scientists worldwide use computers ranging from simple desktop machines to massively parallel arrays of processors to obtain solutions of these equations for the purpose of investigating electromagnetic wave guiding, radiation, and scattering phenomena and technologies. In this chapter, we discuss prospects for using numerical solutions of Maxwell’s equations, in particular the finite-difference time-domain (FDTD) method.

There are seven primary reasons for the expansion of interest in FDTD and related computational solution approaches for Maxwell’s equations [20]:

1. FDTD use no linear algebra. Being a fully explicit computation, FDTD avoids the difficulties with linear algebra that limit the size of frequency-domain integral-equation and finite-element electromagnetics models to generally fewer than 106 electromagnetic field unknowns. FDTD models with as many as 109 field unknowns have been run; there is no intrinsic upper bound to this number.

2. FDTD is accurate and robust. The sources of error in FDTD

calculations are well understood, and can be bounded to permit accurate models for a very large variety of electromagnetic wave interaction problems.

(44)

3. FDTD treats impulsive behavior naturally. Being a time-domain

technique, FDTD directly calculates the impulse response of an electromagnetic system. Therefore, a single FDTD simulation can provide either ultra wideband temporal waveforms or the sinusoidal steady-state response at any frequency with the excitation spectrum.

4. FDTD treats nonlinear behavior naturally. Being a time-domain technique, FDTD directly calculates the nonlinear response of an electromagnetic system.

5. FDTD is a systematic approach. With FDTD, specifying a new structure to be modeled is reduced to a problem of mesh generation rather than the potentially complex reformulation of an integral equation. For example, FDTD requires no calculation of structure-dependent Green’s functions.

6. Computer memory capacities are increasing rapidly. While this trend positively influences all numerical techniques, it is of particular advantage to FDTD methods which are founded on discretizing space over a volume, and therefore inherently require a large random access memory.

7. Computer visualization capabilities are increasing rapidly. While this trend positively influences all numerical techniques, it is of particular advantage to FDTD methods which generate time-marched arrays of field quantities suitable for use in color videos to illustrate the field dynamics.

In fact, the large-scale solution of Maxwell’s equations for electromagnetic wave phenomena using FDTD and similar grid-based time-domain approaches may be fundamental to the advancement of the

(45)

ultra complex and the ultra fast. Maxwell’s equations provide the physics of electromagnetic wave phenomena from do to light, and their accurate modeling is essential to understand high-speed signal effects having wave-transport behavior. A key goal is the computational unification of electromagnetic waves; charge transport in transistors, Josephson junctions, and electro-optic devices; surface and volumetric wave dispersions; and nonlinearities due to quantum effects. Then these can attack a broad spectrum of important problems to advance electrical and computer engineering and directly benefit our society.

2.3.1 Finite-Difference Expressions for Maxwell’s Equations

Imagine a region of space where which contains no flowing currents or isolated charges. Maxwell's curl equations in can be written in Cartesian coordinates as six simple scalar equations. Two examples are:

(2.3.1)

1

y x z y x z

E

H

E

t

z

y

H

E

H

t

y

µ

ε

∂

⎛

⎞

∂

₌

₋

∂

⎜

⎟

∂

_⎝

∂

_⎠

∂

z

∂

=

_⎜

−

_⎟

∂

_⎝

∂

_⎠

⎛

⎞

∂

The other four are symmetric equivalents of the above and are obtained by cyclically exchanging the x, y, and z subscripts and derivatives.

Maxwell’s equations describe a situation in which the temporal change in the E field is dependent upon the spatial variation of the H field,

and vice versa. The FDTD method solves Maxwell's equations by first discretizing the equations via central differences in time and space and then numerically solving these equations in software.

(46)

The most common method to solve these equations is based on Yee's mesh [21] and computes the E and H field components at points on a grid

with grid points spaced ∆x, ∆y, and ∆z apart. The E and the H field

components are then interlaced in all three spatial dimensions as shown in Fig. 2-9. Furthermore, time is broken up into discrete steps of ∆t. The E

field components are then computed at times t = n ∆t and the H fields at

times t = (n+1/2) ∆t, where n is an integer representing the compute step. For example, the E field at a time t = n ∆t is equal to the E field at t =

(n-1) ∆t plus an additional term computed from the spatial variation, or curl, of the H field at time t.

(47)

This method results in six equations that can be used to compute the field at a given mesh point, denoted by integers i, j, k. For example, two of the six are: 1/ 2 1/ 2 ( , 1/ 2, 1/ 2) ( , 1/ 2, 1/ 2) ( , 1/ 2, 1) ( , 1/ 2, ) ( , 1, 1/ 2) ( , , 1/ 2) 1 ( 1/ 2, 1, 1) ( 1/ 2, 1, 1) 1/ 2 ( 1/ 2, 1/ 2, 1) ( 1

(

)

(

)

(

n n x i j k x i j k n n y i j k y i j k n n z i j k z i j k n n x i j k x i j k n z i j k z i

H

t

E

z

t

E

y

E

t

H

y

µ

ε

+ − + + + + + + + + + + + + + + + + + + + + + +

=

∆

+

−

∆

−

∆

=

∆

+

−

∆

1/ 2 / 2, 1/ 2, 1) 1/ 2 1/ 2 ( 1/ 2, , 1/ 2) ( 1/ 2, , 3/ 2)

)

(

)

n j k n n y i j k y i j k

t

H

z

ε

+ − + + + + + + +

∆

−

∆

(2.3.2)

These equations are iteratively solved in a leapfrog manner, alternating between computing the E and H fields at subsequent ∆t/2 intervals.

2.3.2 Stability in FDTD algorithm

Now we discuss that how we determine the time step ∆t. An EM wave propagating in free space cannot go faster than the speed of light. To propagate of one cell requires a minimum time of . When we get to two-dimensional simulation, we have to allow for the propagation in the diagonal direction, which brings the time requirement to 0

/

t

x c

∆ = ∆

0

/( 2 )

t

x

c

∆ = ∆

_{. Furthermore,}

∆ = ∆

t

x

/( 3 )

c

₀ is required in

(48)

three-dimensional simulation. This is summarized by the well-known “Courant Condition” [22] 0

金屬孔洞陣列造成增強性兆赫輻射穿透時孔洞內介質的效應

國 立 交 通 大 學

光 電 工 程 研 究 所

碩 士 論 文

金屬孔洞陣列造成增強性兆赫輻射穿透時

孔洞內介質的效應

Effects of Hole material on Enhanced Terahertz

Transmission through Metallic Hole Arrays

研 究 生：羅 誠

指 導 老 師：潘 犀 靈 教授