二維完美匹配層應用於多重解析法之光學微影模擬

(1)

國立交通大學

電子工程學系電子研究所碩士班

碩士論文

二維完美匹配層應用於多重解析法

之光學微影模擬

Two-Dimensional Perfectly Matched

Layer Applied to Multiresolution

Time-Domain Method for Optical

lithography simulation

研究生 ﹕曾建彰

指導教授 ﹕羅正忠博士

(2)

二維完美匹配層應用於多重解析法

之光學微影模擬

Two-Dimensional Perfectly Matched

Layer Applied to Multiresolution

Time-Domain Method for Optical

lithography simulation

研究生 ﹕曾建彰 Student ﹕Jian-Zhang Zeng

指導教授 ﹕羅正忠博士 Advistor ﹕Dr. Jen-Chung Lou

國立交通大學

電子工程學系電子研究所碩士班

碩士論文

A Thesis

Submitted to Institute of Electronics

College of Electrical and Computer Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of Master of Science

in Electronic Engineering

July 2007

Hsinchu, Taiwan, Republic of China

(3)

二維完美匹配層應用於多重解析法

之光學微影模擬

研究生:曾建彰指導教授:羅正忠

國立交通大學

電子工程學系電子研究所碩士班

摘要

隨著光微影技術向極限推進時，極化效應正變得更加明顯[20]。由於橫向電波和橫向磁波之間結合的不同造成在圖罩(photomask)表面的散射及在晶圓的極化使得成像品質降低[16]。因此我們需要以馬克斯威爾方程式快速且精確的計算晶圓裡的成像。著眼此點，我們企圖將完美匹配層(PML)應用在多重解析法。於是我們先找出一對 Battle-Lemarie 系的尺度和小波函數。在文獻中[1][2]，使用這個可以當作完整基底函數的集合被稱為多重解析度分析。因此我們表明，將多重解析法應用於離散馬克斯威爾方程式的新計畫稱為 MRTD。在此篇論文中，為了簡化我們將只使用尺度函數來開發模擬程式(S-MRTD)。除此之外我們將使用此程式搭配完美吸收邊界層(PML)來模擬開放的二維相移圖罩結構的繞射電磁分析。

(4)

Two-Dimensional Perfectly Matched Layer

Applied to Multiresolution Time-Domain

Method(MRTD) for optical lithography

simulation

Student: Jian-Zhang Zeng Advisor: Dr. Jen-Chung Lou

Department of Electronics Engineering and Institute of Electronics

National Chiao Tung University, Hsinchu, Taiwan

Abstract

Polarization effects are becoming more pronounced [20] as optical lithography technique is pushed towards its limit. Photomask topography scattering, and wafer polarization effects caused by differences in coupling between transverse electric and transverse magnetic waves degrade image quality.[16] Therefore we need to fast and rigorous calculation Latent images in wafer by Maxwell’s equations.

To solve this problem,we attempt to apply’’ Perectly Matched Layer’’ (PML) to ‘’Multiresolution Time-Domain Method’’. Thus we find a pair of cubic spline Battle-Lemarie scaling and wavelet functions. In literature[1][2],the use of these functions as a complete set of basis functions is called multiresolution analysis. Thus we show that the application of multiresolution analysis in the method of moments for the discretization of Maxwell’s equations leads to new multiresolution time-domain

(5)

(MRTD) schemes. In this paper

,

for simplicity

,

we will use the scaling function scheme(S-MRTD) only to develop simulation program. In addition, this program is applied to the rigorous simulation of diffraction from two-dimension open

phase-shifting mask structure. Open boundaries are simulated by the use of a novel formulation of the perfect matching layer(PML) absorber.

(6)

誌

謝

兩年的碩士學業很快就要劃下句號

兩年的碩士學業很快就要劃下句號了

了

了，

了

，

，首先我想感謝我的

首先我想感謝我的

指導教授羅正忠博士

指導教授羅正忠博士，

，

，由於您的細心指導與教誨

，

由於您的細心指導與教誨

由於您的細心指導與教誨，

由於您的細心指導與教誨

，

，讓我在研究

，

讓我在研究

上與許多做人處事上都有莫大的收穫

上與許多做人處事上都有莫大的收穫，

，

，真的使我受益匪淺

真的使我受益匪淺

真的使我受益匪淺，

，

，在

，

在

這裡對老師致上我內心最高的謝意

這裡對老師致上我內心最高的謝意，

，

，感謝您兩年來的教導

，

感謝您兩年來的教導

感謝您兩年來的教導。

。

此外

此外，

，

，我還要感謝張仲興與胡國信學長的指導以及

我還要感謝張仲興與胡國信學長的指導以及

我還要感謝張仲興與胡國信學長的指導以及您所留

您所留

下來的寶貴資源讓我在研究上得以順利解決許多難題及模擬實

驗

驗，

，

，謹此致上最誠摯的謝意

謹此致上最誠摯的謝意

謹此致上最誠摯的謝意。

。

。同時也要感謝曾經一起同甘苦的

同時也要感謝曾經一起同甘苦的

同學們

同學們: 人

人

人傑

傑

傑、

傑

、

、忠

忠樂

忠

樂

樂、

、

、彥銘

、

彥銘

彥銘、

彥銘

、

、德安

、

德安、

德安

、

、智仁

、

智仁

智仁、

、

、睿龍

睿龍

睿龍、

睿龍

、宏仁

、

宏仁

宏仁、

宏仁

、

、俐婷

、

俐婷

俐婷、

、

大峰

大峰、

、

、信智

信智

信智、

、

、勝凱

、

勝凱，

勝凱

，

，謝謝你們在生活及精神上給予的支持

謝謝你們在生活及精神上給予的支持

謝謝你們在生活及精神上給予的支持。

謝謝你們在生活及精神上給予的支持

。

。也

。

也

要感謝跟我住在一起的室友們

要感謝跟我住在一起的室友們:士銘

士銘

士銘、

、

、小誠

小誠

小誠、

、

、飛龍師

飛龍師、

飛龍師

、

、小

小

小 P

PP

P，

，

，謝

謝

謝你們豐富我的生活

謝你們豐富我的生活。

。

。另外也要感謝彤彤姑娘的支持與鼓勵讓

。

另外也要感謝彤彤姑娘的支持與鼓勵讓

我有勇氣克服萬難

我有勇氣克服萬難。

。

最後我要感謝父母苦心的栽培與一再的鼓勵

最後我要感謝父母苦心的栽培與一再的鼓勵，

，

，以及你們無怨

，

以及你們無怨

無悔的付出讓我無後顧之憂完成我的碩士學業

無悔的付出讓我無後顧之憂完成我的碩士學業。

。

。在此獻上我內

。

在此獻上我內

心最深的感恩

心最深的感恩，

，

，辛苦你們了

，

辛苦你們了

辛苦你們了：

：

：爸

爸

爸、

、媽

、

媽

媽、

媽

、

、以及姐姐

以及姐姐

以及姐姐。

。

(7)

Chapter 1 Introduction

1.1 Introduction to Optical Lithography……….……….……1

1.2 Motivation………..……….……….…..3

1.2.1Overview of Optical Imaging System Modeling…………..…....3

1.2.2 Rigorous Electromagnetic Simulation……….…....4

1.2.3 Absorbing Boundary Conditions………...……..5

1.2.4 Numerical Solution of Maxwell’s Equations by MRTD……….5

1.3 Organization of This Thesis……….….….5

Chapter 2 Boundary Conditions for the Finite-Difference

Time-Domain(FDTD) Method

2.1 The FDTD method and the Yee equations……….………..….11

2.2 A Perfectly Matched Layer for the absorption of Electromagnetic Waves………14

2.2.1-The theory and derivation of the PML medium………..…...14

(a) Definition of the PML medium……….…..14

(b) Propagation of a plane wave in a PML medium………...15

(c) Transmission of a Wave through PML-PML Interfaces………...18

2.3 Numerical Stability Condition……….……….………..23

(8)

Chapter 3 Using Battle-Lemarie Scaling Function Based

Multiresolution time-domain(MRTD) schemes

3.1 Fundamentals Of Multiresolution Analysis…….……...29

3.2 Derivation of the S-MRTD Scheme………...30

3.3 Numerical Stability condition………....35

3.4 Numerical Dispersion……….36

Chapter 4 Simulation and Conclusion

4.1 S-MRTD formulation………42 4.2 Simulation results……….……….44 4.3 Conclusion……….………45

References…………

……….52

(9)

Figure Captions

Chapter 1

Figure.1-1 Integrated circuit fabrication process………...7 Figure.1-2 Frequently-encountered circuit geometries……….……….7 Figure.1-3 (a) Optical diffraction blurs dark-bright transitions in the layout. (b)

Photoresist nonlinearity can turn the rather sloppy image into vertical profiles………7 Figure.1-4 Schematic of an exposure system………8 Figure.1-5 High-NA vs low-NA………8 Figure.1-6 Aerial images of a 1:1 line-space pattern with a period of

0.7 / NA

λ

. The

NA is 0.965……….9 Figure.1-7 Latent image contrast loss is less severe than aerial image

degradation………..9 Figure.1-8 Immersion lithography rekindles the need to examine wafer polarization

effects………..……….10 Figure.1-9 A simplified Model for an Imaging System………….…...………….…10

Chapter 2

Figure.2-1 Maxwell’s equations are solved over a cubic grid using the FDTD method. The field components are staggered over the grid ………25 Figure.2-2 (a)The electric field

E

_z(i,j,k) is calculated by summing up the

magnetic field values of the four neighboring nodes. The magnetic field components are assumed to be constant along the line segments 1-2,2-3,3-4,and4-1,and the electric field

E

_z is assumed to be constant over the square surface bounded by 1-2-3-4

(b) The magnetic field Hz(i+1/2,j+1/2,k+1/2) is calculated by summing up the electric field values of the four neighboring nodes. The electric field components are assumed to be constant along the line segments 1-2,2-3,3-4,and4-1,and the magnetic field Hz is assumed to be constant over the square surface bounded by 1-2-3-4.

(10)

Figure.2-3 The transverse electric problem……….….……….26

Figure.2-4 Interface lying between two PML media………….………..……..26

Figure.2-5 The PML technique………..….….…..27

Figure.2-6 Upper-right part of the FDTD grid………...….……..27

Figure.2-7 graphs results obtained using this procedure that illustrate the variation of the numerical phase velocity with propagation angle in a two-dimensional FDTD grid………..……….28

Chapter 3

Figure.3-1 Cubic spline Battle-Lemarie scaling function in space domain……..38

Figure.3-2 Cubic spline Battle-Lemarie scaling function in spectral domain…...38

Figure.3-3 Cubic spline Battle-Lemarie wavelet function in space domain…….39

Figure.3-4 Cubic spline Battle-Lemarie wavelet function in spectral domain….39 Table I Coefficients a(i)……….……40

Figure.3-5 Phase error of MRTD employing scaling functions only: Graphs the variation of the numerical phase error in degrees per wavelength as a function of the parameter, 1/s(1/C), for the case

n =

_a

12

……….………41

Chapter 4

Figure.4-1 The PML regions are labeled as shown in the diagram……….…..46

Figure.4-2 MRTD simulation result……….….…46

Figure.4-3 MRTD simulation result……….……….…47

Figure.4-4 MRTD simulation result……….….47

Figure.4-5 MRTD simulation result……….…48

(11)

Figure.4-7 A new contact hole 2D-PSM design………..….49 Figure.4-8 E-D curve of 100nm contact-hole by the 2D-PSM design………….…49 Figure.4-9 A new contact hole 3D-PSM design………...50 Figure.4-10 E-D curve of 100nm contact-hole by the 3D-PSM design……….…..50 Table 1: FDTD(left) versus MRTD(right) for 2D Quartz object……….…….51 Figure.4-11 Dispersion characteristics of FDTD and S-MRTD, Q=

c t

∆ ∆

/

x

.….51

(12)

Chapter 1 Introduction

1.1 Introduction to Optical Lithography

In the manufacture process,[3] the different layout layers are delineated sequentially onto the silicon wafer by photolithography[4]. This module is the predominant enabler for the 40% per year increase in transistor count per chip. Illustrated in Figure.1-1, the photolithography process starts by coating the silicon substrate with a layer of photoresist. Transistors or wiring patterns are then imaged onto the photoresist by optical lithography. After development of the resist latent image, materials are removed from or deposited onto the wafer by processing steps such as etching, deposition, chemical mechanical polishing (CMP), and ion

implantation. The photoresist is subsequently removed, leaving the substrate ready for another cycle of processing steps. Fabrication is complete after cycling a silicon wafer through these processing modules 30 to 40 times.

Circuit geometries printed by optical lithography can be generalized into the two classes illustrated in Figure.1-2(a): line-space and contacts. As a simplification, lines and spaces are usually assumed periodic with period p and dimension d

[Figure.1-2(b)]. For example, image of a periodic line-space array with a 1 : 1 duty ratio is shown in Figure.1-3(a). Ideally transition from dark to bright regions should be sharp so that the photoresist is either exposed or unexposed. But optical diffraction blurs the transition. Fortunately photoresist development nonlinearity can turn the rather sloppy images into vertical profiles, as illustrated in Figure.1-3(b).

(13)

The photolithography imaging system is shown schematically in Figure.1-4. A narrowband light source of wavelength

λ

is placed at the focal plane of the condenser represented by lens L1, illuminating the photomask by Ko”hler's

method[5]. Image of the photomask, template that contains the layout geometries, is formed by the projection optics onto the wafer. The projection optics is represented by the two lenses L2 and L3 in the Figure1-4; in actual systems it can consist of more than 30 lenses[6]. Because of the Fourier transforming property of lens L2[7], energy transmitted through the photomask forms a distribution in the pupil plane that is proportional to the mask spectrum. Low-spatial-frequency components pass closer to the center of the pupil whereas higher frequency components are nearer the peripheral of the pupil. The highest frequencies are cut off by the pupil. Observing from the wafer plane, the image-forming rays are restricted in angular extent by the pupil, subtending an angle

θ

at the wafer. The sine of the angle

θ

is the numerical aperture (NA) of the system: NA = sin

θ

. In the pupil plane an image of the source is formed. For a circular source and pupil in conventional optical lithography, the partial coherence factor

σ

is defined as the ratio between the sizes of the source image (the effective source) and the pupil. For the system drawn in Figure1-4.

a

b

σ

=

.

Mathematically the imaging process can be described by[5]

2 [( ' '') ( ' '') ]

( , )

( , ) (

',

')

(

'',

'')

( ', ')

( '', '')

i f f x g g y

'

''

,

I x y

J f g P f

f g

g P

f

g

O f g O

f

g e

π

df dg df dg dfdg

+∞ _∗ −∞ ∗ + − + −

=

_∫

⋅⋅⋅

_∫

+

⋅

where

I x y

( , )

is the image intensity at location (x,y) on the wafer,

J

is the

effective source,

P

is the pupil function,

O

is the mask spectrum, and the asterisk

(14)

Similar to Rayleigh's resolution criterion, resolution of an exposure

system scales with (

λ

/NA). It is customary to represent the critical dimension (CD) of the printed feature size as

1

,

CD

k

NA

λ

= ⋅

where k1 is a process-related factor measuring the ease of lithography. Smaller dimensions can be printed by decreasing the wavelength, increasing the numerical aperture, and reducing k1. All three measures must be applied to sustain Moore's law and to achieve the ultimate resolution.

1.2 Motivation

1.2.1 Overview of Optical Imaging System Modeling

A modern, complex optical imaging system can have several lens, polarizer, mask, aperture,wafer, and pupil components and is usually very large in relation to the wavelength of light. Thus it is not practical to simulate the imaging system with the FDTD method .However, the science of optics can be applied to the imaging system components for modeling purposes. The Hopkins theory of partially coherent imaging is commonly used to calculate aerial images in lithography. .

To print or inspect smaller features, optics designers build optics with higher numerical apertures(Figure.1-5). For numerical apertures exceeding 0.5 to 0.6, the paraxial approximations made in scalar imaging theory are invalid and theory has been extended by Daniel C.Cole et al.[8] who removed the paraxial ray approximation in the projection optic model. This led to the “Radiometric Correction Factor” that extended the usefulness of the projection optic model to numerical apertures in the

(15)

range 0.6 to 0.7. But even with the radiometric correction factor, the scalar theory is deficient for high numerical aperture lithography ( NA) due to the differences in the way TE and TM plane waves behave, especially at the image plane where highly oblique plane waves exist. M S.Yeung[9] realized that a vector formulation is required for high NA imaging inside thin film stacks and generalized the Hopkins’ formula for vector fields. But he then points out that, upon entering the photoresist, the plane waves bend towards the normal reducing their degree of obliquity, and reverts back to a scalar theory which he claims is “sufficient” for the practical simulation of aerial images in planar thin-film structures for numerical apertures at least as high as 0.6. But the recent interest in immersion lithography rekindles the need to examine wafer polarization effects. With the incident medium being a liquid of refractive index

n

liquid higher than 1, resolution of immersion lithography is CD= ₁ liquid

k

n

NA

λ

⋅

For 193-nm lithography, a promising immersion liquid is water having a refractive index of 1.43, giving a 43% resolution improvement. With the use of an immersion liquid, ray angles in resist becomes larger, resulting in more severe polarization effects.(see Figure.1-6,7,8 we can known)

1.2.2 Rigorous Electromagnetic Simulation

As device dimensions are pushed deeper and deeper into the submicron regime ,polarization and other complex electromagnetic effects that arise because of mask topography are becoming increasingly more important to model.

(Figure.1-9) The science of Fourier Optics is applied to model the source, illumination optic and projection optic. Rigorous ElectroMagnetic

(16)

Simulation (REMS) may be required to model the object (photomask). The formation of the aerial image inside a film stack(Latent image) at the image plane may be model either by REMS or by thin-film-stack theory.[11]

1.2.3 Absorbing Boundary Conditions

In 1994, Berenger published a paper[10] on a new type of boundary condition which he named “Perfectly Matched Layers” or PML. The boundary condition is essentially a non-physical material with the special property that it absorbs electromagnetic radiation without refection for all frequencies and angles of

incidence. In this paper, the new boundary condition was implemented into S-MRTD and is explained in Chapter 2.

1.2.4 Numerical Solution of Maxwell’s Equations by MRTD

In this paper, the recently developed multiresolution time-domain (MRTD) method[12] is applied to fast and rigorous mask diffraction simulation, to overcome the limitations of FDTD. An outline of the MRTD formulation is described in Chapter 3. In addition, In chapter 4, using a 2D scattering example to allow the use of an extremely fine computational mesh to achieve good convergence.

1.3 Organization of This Thesis

The thesis includes four chapter. In chapter 1, we make an introduction to describe the background and polarization effects of the optical lithography and the role in the semiconductor industry that optical lithography is playing. In addition, we make a

(17)

brief to describe simple concept of the imaging process mathematically. Then, we describe the motivation of this thesis. Finally, we introduce PML and MRTD techniques respectively.

In chapter2, we will first specify a numerical scheme for solving the Maxwell’s equations, used by Yee.[13] Yee was one of the first to replace Maxwell’s equations by a set of finite difference equations. Then, we will illustrate PML technique for two-dimensional problems.

In chapter3, we first explain the most important properties of the multiresolution analysis. Then, we will only derive the S-MRTD scheme for a homogeneous medium Finally,in chapter4, we will use S-MRTD and PML schemes to rigorous simulate two-dimension phase-shifting mask structure by using matlab. Then, we will make a conclusion with simulation in the whole thesis.

(18)

Figure.1-1 Integrated circuit fabrication process

( )

a

( )

b

Figure.1-2 Frequently-encountered circuit geometries. (a) lines, spaces, and contacts (b) periodic line-space pattern

( )

a

( )

b

Figure.1-3 (a) Optical diffraction blurs dark-bright transitions in the layout. (b) Photoresist nonlinearity can turn the rather sloppy image into vertical profiles.

(19)

Figure.1-4 Schematic of an exposure system.

(20)

Figure.1-6 Aerial images of a 1:1 line-space pattern with a period of

0.7 /

λ

NA

. The NA is 0.965.

(21)

Figure.1-8 Immersion lithography rekindles the need to examine wafer polarization effects.

(22)

Chapter 2 Boundary Conditions for the Finite-Difference

Time-Domain(FDTD) Method

2.1 The FDTD method and the Yee equations

The continuous form of the Maxwell equations for linear, isotropic, non-magnetic, nondispersive materials are written:

D H J t ∂ ∇ × = + ∂

(1.a) B E t ∂ ∇ × = − ∂ (1.b) with the constitutive relations:

D

E

B

H

J

E

ε

µ

σ

=

(2)

Where

ε µ

, ,

and

σ

are, respectively, the permittivity, permeability, and

conductivity of the material. In general, the parameters

ε µ

,

, and

σ

are functions of the frequency of the electromagnetic wave. For the application in hand, however,

they are assumed to be constant because of monochromatic excitation. Using Stokes' theorem, (1.ab) and (2) can be re-written in the weak form

}

l s

D

H dl

J dS

t

∂

=

_

+

∂



∫

i

(3) l s

B

E dl

dS

t

∂

= −

∂

∫

i

(4)

(23)

Where

l

F dl

∫

i

. and

s

F dS

∫

i

_{represent, respectively, the line integral and}

surface integral of a variable

F

. (1.ab) and (2) also are equivalent to the following system of scalar equations:

,

y x z

E

B

E

t

y

z

∂

−

=

−

∂

,

y x z

B

_E

t

z

x

∂

−

=

−

∂

,

y z x

E

B

E

t

y

x

∂

=

−

∂

,

y x z x

H

D

H

J

t

y

z

∂

=

−

∂

,

y x z y

D

H

J

t

z

x

∂

=

−

∂

.

y z x z

H

D

H

J

t

x

y

∂

=

−

∂

(5)

Following the FDTD method proposed by Yee [6] in which the differential

δ

t

is replaced by

∆

t

and by

∆

x

, (3)(4) or (5) are solved using a cubic grid in which the field components are staggered and occupy distinct locations in space as shown in Fig. 2.1 The surface integral and line integral are thus evaluated on square surfaces as shown in Fig. 2.2. This discretization scheme[17] leads to three scalar equations in two-dimensional analysis for the transverse electric (TE) or the transverse magnetic (TM) polarization, but six equations in three-dimensional analysis for the field components

E E E H H

x

,

y

,

z

,

x

,

y

,

and

H

z

:

(24)

Where

α

=

(2

ε σ

−

t

) /(2

ε σ

+

t

),

β

=

(

t

/

x

) [2 /(2

⋅

ε σ

+

t

)],

and

(

t u x

/

)

γ

=

. The superscript of the field variables stands for the time step (time =

n

⋅ ∆

t

), the subscript represents the direction of the field, and (i, j , k)

(25)

signifies the node position at

(

i x j y k z

∆

,

∆

,

∆

)

.

2.2 A Perfectly Matched Layer for the

absorption of Electromagnetic Waves

2.2.1 The theory and derivation of the PML medium

(a) Definition of the PML medium

We will set the equation of a PML medium for two-dimensional problems,first in the TE(transverse electric) case. In Cartesian coordinates let us consider a problem that is without variation along z,with the electric field lying in the (x,y) plane(Fig.2.3) . The electromagnetic field involves three components only,Ex,Ey,Hz,and the Maxwell equations reduce to a set of three equations. In the most general case,which is a medium with an electric conductivity

σ

and a magnetic conductivity

σ

∗

,these equations can be written as

0 0 0

*

x z x y z y y z x z

E

H

E

t

y

E

H

E

t

x

E

H

E

u

H

t

y

x

ε

σ

ε

σ

∂

+

=

∂

_∂

+

= −

∂

+

=

−

∂

(1)

Moreover,if the condition

0 0

*

u

σ

(26)

is satisfied,then the impendence of the medium (1) equals that of vacuum and no reflection occurs when a plan wave propagates normally across a vacuum-medium interface.

We will now define the PML medium in the TE case. The cornerstone of this definition is the break of the magnetic component Hz into two subcomponents which we will denote as

H

_zxand

H

_zy. In the TE case,a PML medium is defined as a medium in which the electromagnetic field has four

components,

E E H

_x

,

_y

,

_zx

,

H

_zy,connected through the four following equations:

0 0 0 0

(

)

(

)

*

z x zy x y x y z x zy x y y zx x zx zy _x y zy

H

E

t

y

E

H

E

t

x

E

H

u

H

t

x

H

E

u

H

t

y

ε

σ

ε

σ

∂

+

∂

+

=

∂

+

= −

∂

+

= −

∂

_∂

+

=

∂

(3)

where the parameters

(

σ σ

_x

,

_x

*,

σ σ

_y

,

_y

*)

are homogeneous to electric and magnetic conductivities.

(b) Propagation of a plane wave in a PML medium

Let us consider a wave whose electric field of magnitude

E

₀ forms an angle

ϕ

with the y axis (Fig2.3). We will denote as

H

_zx₀and

H

_zy₀ the magnitudes of the magnetic components

H

_zxand

H

_zy. If a plane wave propagates in the PML

(27)

( ) 0

sin

iw t x y x

E

= −

E

ϕ

e

−α −β ( ) 0

cos

iw t x y y

E

=

E

ϕ

e

−α −β ( ) 0 iw t x y zx zx

H

=

H

e

−α −β ( ) 0 iw t x y zy zy

H

=

H

e

−α −β

(4) Where w is the pulsation of the wave, t is the time,and

α

and

β

are complex constants. Since the magnitude

E

₀ is given,the set of Eq.(4)involves four unknown quantities to be determined,

α β

, ,

H

_zxo

,

H

_zyo. Enforcing

E E H

_x

,

_y

,

_zx

,

H

_zy from(4) in the PML equations (3) yields the following set of equations connecting the four unknowns: 0 0

sin

0

sin

(

0 0

)

y zx zy

E

i

E

H

w

σ

ε

ϕ

−

ϕ β

=

+

(5) 0 0

cos

0

s cos

(

0 0

)

x zx zy

E

i

E

H

w

σ

ε

ϕ

−

ϕ α

=

+

(6) 0 0 0 0

*

cos

x zx zx

u H

i

H

E

w

σ

_α

_ϕ

−

=

(7) 0 0 0 0

*

sin

y zy zy

u H

i

H

E

w

σ

β

ϕ

−

=

(8) Obtaining

H

_zx₀and

H

_zy₀ from (7) and (8) and bringing them respectively into (5) and (6) yields 0 0 0 0 0

cos

sin

(1

)sin

[

]

(1

(

* /

))

1 (

* /

)

y x y

u

i

w

i

u w

i

u w

σ

α

ϕ

β

ϕ

ε

ϕ β

ε

σ

−

=

+

−

(9)

(28)

0 0 0 0 0

cos

sin

(1

) cos

[

]

(1

(

* /

))

1 (

* /

)

x x y

u

i

w

i

u w

i

u w

σ

α

ϕ

β

ϕ

ε

ϕ α

ε

σ

−

=

+

−

(10)

This system of two equations connects the unknowns

α

and

β

. It may be solved by means of writing the ratio (9) over (10),

0 0

1 (

/

)

sin

cos 1

(

/

)

y x

i

w

i

w

σ ε

β

ϕ

α

ϕ

σ ε

−

=

−

(11)

And then obtaining

α

2 from (11) and (10) and

β

2 from (11) and (9). That yields two sets of (

α

,

β

) of opposite signs for two opposite directions of propagation. Choosing the positive sign we have

0 0 0

(1

x

) cos

u

i

G

w

ε

σ

α

ϕ

ε

=

−

0 0 0

(1

y

)sin

u

i

G

w

σ

ε

β

ϕ

ε

=

−

(12) Where G=

w

_x

cos

2

ϕ

+

w

_y

sin

2

ϕ

(13)

0 0

1 (

/

)

1 (

* /

)

x x x

i

w

i

u w

σ ε

σ

−

=

−

(14) 0 0

1 (

/

)

1 (

* /

)

y y y

i

w

i

u w

σ ε

σ

−

=

−

(15)

Denoting as

ψ

any component of the field ,

ψ

₀ its magnitude,and c the speed of light , with (4)and (12) , we can write

0 0 ( sin / ) ( cos / ) ( ( cos sin ) / ) 0

.

y x cG x cG y iw t x y cG

e

ϕ ϕ

e

σ ϕ ε

e

σ ϕ ε

ψ ψ

₌

− + − − (16) The last two unknowns

H

_zx₀ and

H

_zy₀ can be found as functions of

α

and

β

from (7) and (8) , and then enforcing the

α

and

β

values (12) yields

(29)

2 0 0 0 0

1 (

/

)

cos

zx x

H

E

u

w

G

ε

ϕ

=

2 0 0 0 0

1 (

/

)

sin

zy y

H

E

u

w

G

ε

ϕ

=

(17)

Taking into account (13), the summation of

H

_zx₀ and

H

_zy₀ is then

0 0

(

0

/

0

)

H

=

E

ε

u G

(18) and the ratio Z of the electric magnitude over the magnetic one is

0 0

1 /

Z

u

G

ε

=

(19)

For formulas (16) and (19), an important occurrence is when both (

σ σ

_x

,

_x

*

) and (

σ σ

_y

,

_y

*

) satisfy the condition(2). Then, the quantities

w w G ,

_x

,

_y

,

equal unity at any frequency, and so the expression of the wave components (16) and of the impedance (19) become respectively

0 0 ( sin / ) ( cos / ) ( ( cos sin ) / ) 0 y x c x c y iw t x y c e ϕ ϕ e σ ϕ ε e σ ϕ ε

ψ ψ

₌ − − _× − −

(20) 0/ 0. Z = u

ε

(21)

(c) Transmission of a Wave through PML-PML Interfaces

In this section,we will address the problem of a plane wave moving from a PML medium to another one. We will prove that in particular

cases, with

adequate sets of parameters

σ σ

_x

,

_x

*,

σ σ

_y

,

_y

*,

the transmission is perfect and reflectionless at any incidence angle. These particular cases will be the basis of the perfectly matched

(30)

layer.

We first consider the case of two PML media separated by an interface normal to the x axis(fig2-4) . Let us denote

θ

₁ and

θ

₂ the angles of the incident and

transmitted electric fields

E

_i and

E

_t,with respect to the interface plane. In the case of (fig.2) we assume the interface to be infinite and incident wave to be plane.

First,both the reflected and transmitted waves must be plane too. Second, the ratios of these waves over the incident one must be without variation when moving on the interface. So, for any component

ψ

of the incident and transmitted fields, and for two points A and B of the interface, we can write

( )

t t i i

B

A

B

A

ψ

=

ψ

(22) Denoting as d the distance from A to B, and

G

₁ ,

G

₂ as the quantities(13) of each media, with (16) we have

1 1 1 1 0 1 ( sin / ) ( sin / )

( )

iw d cG y cG d i

B

i

A e

θ σ θ ε

ψ

− −

=

(23) 2 2 2 2 0 2 ( sin / ) ( sin / )

( )

iw d cG y cG d t

B

t

A e

θ σ θ ε

ψ

− −

=

(24)

Since (22) is true for any distance d, the exponential factors of (23) and (24) must be equal. So the relation

1 ₁ 2 ₂ 0 1 0 2

sin

(1

i

y

)

(1

i

y

)

w

G

w

G

σ

θ

σ

θ

ε

−

= −

(25)

Must be satisfied, where

2 2

cos

sin

k xk k yk k

(31)

Let us now consider the incident, reflected, and transmitted electric and magnetic fields

E E E H H H

_i

,

_r

,

_t

,

_i

,

_r

,

_t(fig.2.4). First, continuity of the fields

E

_y and

zx zy

H

+

H

lying in the interface yields the following set of equations:

1 1 2

cos

i r t

E

θ

−

E

θ

=

E

θ

(27)

i r t

H

+

H

=

H

(28) Second, denoting as

E E

_io

,

_ro

,

E

_to

,

the magnitudes of

E E E

_i

,

_r

,

_t

,

and setting x=0 in the interface, with (2.2.11) and (2.2.14) we can write

1 1 1 ( sin / )(1 ( y / o )) iw y cG i w iwt i io

E

=

E e

− θ − σ ε

e

(29) 1 1 1 ( sin / )(1 ( y / o )) iw y cG i w iwt r ro

E

=

E e

− θ − σ ε

e

(30) 2 2 2 ( sin / )(1 ( _y / _o )) iw y cG i w iwt t to

E

=

E e

− θ − σ ε

e

(31) 1

/

,

i i

H

=

E Z

H

_r

=

E

_r

/

Z

₁

,

H

_t

=

E Z

_t

/

₂

.

(32)

As a consequence of the Snell-Descartes law (25),(26),the three exponentials on space of (29),(30),(31) are equal. So, reporting

E

_i….,

H

_t from (29),(30),(31) into (27),(28) this system (27),(28)becomes

1 1 2

cos

io ro to

E

θ

−

E

θ

=

E

θ

(33) 1 1 2 io ro to

E

Z

+

Z

=

Z

(34)

Defining the reflection factor as the ratio of the electric components in the interface, that is,

−

E

_ro

cos

θ

₁

/

E

_io

cos

θ

₁ ,and then solving the set (33),(34) for this ratio, the reflection factor

r

_p for the TE case is

(32)

2 2 1 1 2 2 1 1

cos

p

Z

r

Z

θ

−

=

+

(35)

With (19) , we can rewrite (35) as

1 2 2 1 1 2 2 1

cos

p

G

r

G

θ

−

=

+

(36)

where

G

₁ and

G

₂ are functions of

θ

₁ and

θ

₂ through (26)

Let us consider an interface lying between media having same

σ

_y and

σ

_y∗ conductivities, that is, a

(

σ σ σ σ

_x₁

,

_x∗₁

,

_y

,

_y∗

)

and

(

σ σ σ σ

_x₂

,

_x∗₂

,

_y

,

_y∗

)

media. Then the (25) becomes 1 2 1 2

sin

G

θ

=

(37)

If, moreover, the two media are matched ones, that is,

(

σ σ

_x₁

,

_x∗₁

),(

σ σ

_x₂

,

_x∗₂

)

, and

(

σ σ

_y

,

∗_y

)

satisfying (2), we have

G

₁

=

G

₂

=

1

, so (37) reduces to

1 2

θ θ

=

(38) And the reflection factor (36) is then

0

p

r =

(39) The general frame of the PML technique is pointed out on Fig.2.5.

As a summary of this section we can now say that the reflection factor between two PML media whose conductivities satisfy (2) is null: These are

the reflection factor is null at an interface normal to x lying between a vacuum and a

(

σ σ

_x

,

_x∗

,0,0)

matched medium or between a

(0,0,

σ σ

_y

,

_y∗

)

and a

(33)

the reflection factor is null at an interface normal to y lying between a vacuum and a

(0,0,

σ σ

_y

,

∗_y

)

matched medium or between a

(

σ σ

_x

,

_x∗

,0,0)

and a

(

σ σ σ σ

_x

,

_x∗

,

_y

,

_y∗

)

matched media.

Let us consider, for instance, the upper-right part of a gridded computational domain (Fig2-6). With the usual FDTD notations, (3.b) and (3.c) yield the following equations that can be applied in the whole layer, except in the interface for

E

_y(see below), ( ) / ( ) / 1 1/ 2 1/ 2 1/ 2 1/ 2

(1

)

( ,

1/ 2)

( ,

1/ 2)

( )

[

(

1/ 2,

1/ 2)

(

1/ 2,

1/ 2)

(

1/ 2,

1/ 2)

(

1/ 2,

1/ 2)]

x o x o i t i t n n y y x n n zx zy n n zx zy

e

E

i j

e

E i j

i

x

H

i

j

H

i

j

H

i

j

H

i

j

σ ε σ ε

σ

− ∆ − ∆ + + + + +

−

+

=

+

−

×

+

−

+

−

+

(40) ( 1/ 2) / 1 1/ 2 ( 1/ 2) /

(

1/ 2,

1/ 2)

(

1/ 2,

1/ 2)

(1

)

(

1/ 2)

[

(

1,

1/ 2)

( ,

1/ 2)],

x o x o i t u n n zx zx i t u x n n y y

H

i

j

e

H

i

j

e

i

x

E i

j

E i j

σ σ

σ

∗ ∗ − + ∆ + − − + ∆ ∗

+

=

+

−

+

×

+

−

+

(41)

Where

σ

_x and

σ

_x∗ are functions of x(I).(on Fig2-6) For the

E

_y component lying in the interface, the magnetic field has one component

H

_z on one side and two subcomponents

H

_zx

,

H

_zy on the other. That is a result of using the Maxwell equations in the inner volume, but it has no physical significance. The finite

difference equations have to be modified. So, in the right side interface normal to x, (40) becomes

(34)

( ) / 1 ( ) / 1/ 2 1/ 2 1/ 2

( ,

1/ 2)

( ,

1/ 2)

(1

)

( )

[

(

1/ 2,

1/ 2)

(

1/ 2,

1/ 2)

(

1/ 2,

1/ 2)].

x o x o il t n n y y il t x n n n zx zy z

E

il j

e

E il j

e

il

x

H

il

j

E

il

j

H

il

j

σ ε σ ε

σ

− ∆ + − ∆ ∗ + + +

+

=

+

−

×

+

−

+

(42)

2.3 Numerical Stability Condition

In the two-dimension case

2 2

1

_o

2 t

c

x

y

∆

∆ ≤

=

+

and

∆ = ∆ = ∆

x

y

2.4 Numerical Dispersion Relation

In the two-dimension TE(TM) case

2

1

1 sin

sin

2

y x

k

y

w t

k

x

c t

x

y





















_{ =}



_{+ }

_



















_

_



_

_





























where

k

_x and

k

_y are, respectively, the x- and y-components of the numerical wavevector, and w is the wave angular frequency. To quantitatively assess the dependence of numerical dispersion upon FD-TD grid discretization, we shall take as the case, assuming for simplicity square unit cells (

x

=

y

=

) and wave

propagation at an angle

α

with respect to the positive

x-axis(

k

_x =

k

cos

α

;

k

_y

=

k

sin

α

). Then numerical dispersion relation () simplifies to

(35)

2

2 2

cos

2

sin

2

2 w t

k

c t

α



⋅





⋅









_{ =}

₊

























_

_

_

_

By applying the following Newton’s method iterative procedure:

( )

2 2 1

sin

sin(2

)

sin(2

)

i i i i i i

Ak

Bk

C

k

A

Ak

B

Bk

+

−

= −

+

cos

2 A

=

⋅

α

,

sin

2 B

=

⋅

α

, 2 2

sin (

)

2 w t

C

c t





= 







So we can obtain the numerical wavevector

k

_final that the numerical phase velocity

p

v

is given by

p final

w

v

k

=

.

See figure.2-7 graphs results obtained using this procedure that illustrate the

variation of the numerical phase velocity with propagation angle in a two-dimensional FDTD grid.

(36)

Fig2.1 Maxwell’s equations are solved over a cubic grid using the FDTD method. The field components are staggered over the grid.

(a)

(b)

Fig2.2-(a)The electric field

E

_z(i,j,k) is calculated by summing up the magnetic field values of the four neighboring nodes. The magnetic field components are assumed to be constant along the line segments 1-2,2-3,3-4,and4-1,and the electric field

E

_z is assumed to be constant over the square surface bounded by 1-2-3-4.

(b) The magnetic field Hz(i+1/2,j+1/2,k+1/2) is calculated by summing up the electric field values of the four neighboring nodes. The electric field components are assumed to be constant along the line segments 1-2,2-3,3-4,and4-1,and the magnetic field Hz is assumed to be constant over the square surface bounded by 1-2-3-4.

(37)

Fig-2.3The transverse electric problem

(38)

Fig-2.5 The PML technique.

(39)

wave propagation angle

α

Figure.2-7 graphs results obtained using this procedure that illustrate the variation of the numerical phase velocity with propagation angle in a two-dimensional FDTD grid.

(40)

Chapter 3 Using Battle-Lemarie Scaling Function Based

Multiresolution time-domain(MRTD) schemes

3.1 Fundamentals Of Multiresolution Analysis

In literature [1], [2], the use of scaling and wavelet functions as a complete set of basis functions is called multiresolution analysis.The most important properties of the scaling and wavelets functions defining a multiresolution analysis are now

summarized:

1. All scaling and wavelet functions are localized both in the frequency and the space/time domains.

2. The Fourier transform of a scaling function(Fig.3-1) is a low-pass

function(Fig.3-2), while the Fourier transform of the corresponding mother wavelet(Fig.3-3) is a bandpass function(Fig.3-4). This implies that, in a multiresolution expansion, scaling functions can be utilized to sample the low-frequency content of the field/signal, while wavelets sample the high-frequency content.

3. The scaling function of one resolution level, m, are orthonormal to the scaling functions of the same resolution level, but not the scaling functions of other resolution levels:

(41)

,

mn mv nv

φ φ

=

δ

(3-1) This implies that, when using scaling functions in an expansion, we cannot mix different resolution levels.

4. The scaling function are orthogonal to all wavelets of any higher-order resolution level:

φ ψ

_mn

,

_uv

=

0 m

≤

u

(3-2) while the wavelet function of any resolution level are orthonormal to

all other wavelets at any resolution level:

,

mn uv mu nv

ψ ψ

=

δ δ

(3-3) The above two properties allow us to mix scaling functions of one resolution level

with wavelets of one or more resolution levels of the same or higher order, m. 5. Any scaling function of m resolution can be expanded in the form of

wavelets of all lower resolution levels: 1

( )

m mn mn ji ji j i

x

d

ψ

x

φ

− ∞ =−∞ =−∞

=

∑ ∑

(3-4)

3-2 Derivation of the S-MRTD Scheme

The S-MRTD and W-MRTD schemes[12] are derived using cubic spline Battle-Lemarie scaling and wavelet functions .

At first, for simplicity, we will only derive the S-MRTD scheme for a

homogeneous medium. The derivation is similar to that of Yee’s FDTD scheme which uses the method of moments[19] with pulse functions as expansion and test functions. For the derivation of S-MRTD, the field components are represented by a series of scaling functions in space and pulse functions in time.

(42)

2-D (TE) S-MRTD Scheme: Maxwell’s vector equation

E

H

t

ε

∂

∇ ×

=

∂

(1) 0

H

E

u

t

∂

∇ × = −

∂

(2)

for a homogeneous medium with the permittivity

ε

and the permeability

u

₀ may be written in the form of three scalar Cartesian equations as

z x

H

E

y

ε

t

∂

=

∂

(3) y z

E

H

x

ε

t

∂

−

=

∂

(4) 0 y x

E

z

E

H

u

y

x

t

∂

−

=

∂

(5)

The electric and magnetic field components incorporated in these equations are expanded as following 1/ 2, 1/ 2 , , ,

( , )

x

( )

x k l m k l m k l m

E r t

E

φ

h t

φ

x

φ

y

+∞ + + =−∞

=

∑

, 1/ 2 1/ 2 , , ,

( , )

y

( ) ( )

( )

y k l m k l m k l m

E r t

E

φ

h t

φ

x

φ

y

+∞ + + =−∞

=

∑

1/ 2 1/ 2, 1/ 2 1/ 2 1/ 2 1/ 2 , , ,

( , )

z

( )

z k l m k l m k l m

H r t

H

φ

h

t

φ

x

φ

y

+∞ + + + + + + =−∞

=

∑

(6)

(43)

space and time indices related to the space and time coordinates via

x

= ∆

l x

,

y

= ∆

m y

and

y

= ∆

k t

, where

∆

x

,

∆

y

and

∆

t

represent the space and time discretization intervals in x-, y- and t-direction. The function

h x

_m

( )

is defined as

( )

(

)

m

x

h x

h

m

x

=

−

∆

(7)

with the rectangular pulse function

1 1/ 2

( )

1/ 2

0 1/ 2

for x

h x

for x

<





=

_

=



_>



(8)

The function

φ

m

( )

x

_{is defined as}

( )

(

)

m

x

m

x

φ

=

φ

−

∆

(9)

where

φ

( )

x

represents the spline Battle-Lemarie scaling function [13], [14] depicted in Fig. 3-1. Assuming the Fourier transformation

( )

x e

j xλ

dx

φ λ

+∞

φ

−∞

=

_∫

(10) and

1 ( )

( )

2

j x

x

e

λ

d

φ

φ λ

λ

π

+∞ ₋ −∞

=

_∫

(11)

Respectively, the closed-form expression of the scaling function in spectral domain is given by

4 2 4 6

sin( )

1

2 ( )

4

2

4

1 sin ( )

sin ( )

2 ₃

₂

₅

₂

₃₁₅

₂

λ

φ λ

_λ

λ









=

_

_

⋅





₋

₊

₋





(12)

(44)

with the low-pass spectral domain characteristics shown in Fig.3-2.

We insert the field expansions in Maxwell’s equations and sample the equations using pulse functions as test functions in time and scaling functions as test functions in space. For the sampling with respect to time, we need the following integrals[10]

' , '

( )

m m m m

h x h

x dx

δ

x

+∞ −∞

=

∆

∫

(13) where

δ

_{m m}_{, '} represents the kronecker symbol

, '

1 '

0 '

m m

for

m

for

m

δ

_{= }



=

≠



(14) and ' 1/ 2 , ' , ' 1

( )

m m m m m m

h

x

h x

dx

x

δ

+∞ + + −∞

∂

=

−

∂

∫

(15) For the sampling with respect to space, we use the orthogonality

relation for the scaling functions [17]

' , '

( )

.

m

x

m

x dx

m m

x

φ

δ

+∞ −∞

=

∆

∫

(16) To calculate the integral corresponding to (15) for scaling functions, we

make use of the closed form expression of the scaling function in spectral domain. According to Galerkin’s method [7], for complex basis functions, one has to choose the complex conjugant of the basis functions as test functions. We then obtain

2 ' 1/ 2 0

( )

1 ( )

m

( )

sin ( '

1/ 2)

m

x

dx

m

d

x

φ

φ λ λ

λ

π

+∞ ∞ + −∞

∂

=

− +

∂

∫

(17)

where

φ λ

( )

is given by (12). This integral may be evaluated numerically resulting in

(45)

' 1/ 2 , '

( )

m

( )

m m i m i

x

dx

a i

x

φ

+∞

δ

+∞ + + −∞ =−∞

∂

=

∂

∑

∫

(18) The coefficients a(i) for

0 ≤ ≤

i

8

are shown in Table I, and the coefficients a(i) for i < 0 are given by the symmetry relation a(-1-i)=-a(i). The Battle-Lemarie scaling function does not have compact but only exponential decaying support and thus, the coefficients a(i) for i > 8 are not zero. However, we found that these coefficients are negligible, affecting the accuracy of the field computation only for very low values of the wave vector. We therefore use the approximation

8 ' 1/ 2 , ' 9

( )

m

( )

m m i m i

x

dx

a i

x

φ

+

δ

+∞ + + −∞ =−

∂

≈

∂

∑

∫

(19) in order to obtain a MRTD scheme useful for practical applications.

For sake of simplicity,let us demonstrate the sampling of Maxwell’s

equations by using only scaling functions as expansion and testing functions. By applying Galerkin’s method,the system of (3),(4),(5) takes the form:

8 1 1/ 2, 1/ 2, 1/ 2 1/ 2, 1/ 2 9

( )

x x z k l m k l m k l m i i

t

E

a i

H

y

φ φ φ

ε

+ + + + + + + + =−

∆ 



=

+

_

_

∆ 

∑



8 1 , 1/ 2 , 1/ 2 1/ 2 1/ 2, 1/ 2 9

( )

y y z k l m k l m k l i m i

t

E

a i

H

x

φ φ φ

ε

+ + + + + + + + =−

∆ 



=

+

_

−

_

∆ 

∑



8 8 1/ 2 1/ 2, 1/ 2 1/ 2 1/ 2, 1/ 2 1/ 2, , 1/ 2 9 9 0

( )

z z x y k l m k l m k l m i k l i m i i

t

H

a i

E

a i

E

u

φ φ + φ + φ + + + − + + + + + + =− =−

∆ 



=

+

_

−

_

∆ 

∑



(20) where l,m, and k are the discrete space and time indices. In fact, the total field at a particular space point for the S-MRTD scheme may be calculated from the field expansions, see (6) , by sampling them with delta test functions in space and time

二維完美匹配層應用於多重解析法之光學微影模擬

國 立 交 通 大 學

電子工程學系 電子研究所碩士班

碩 士 論 文