分析光柵結構的有限差分法與傅立葉模態法

(1)

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

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

National Taiwan University Master Thesis

分析光柵結構的有限差分法與傅立葉模態法

Finite-Difference and Fourier Modal Methods for the Analysis of Gratings Structures

沈祺凱 Chi-Kai Shen

指導教授：邱奕鵬博士 Advisor: Yih-Peng Chiou, Ph.D.

中華民國九十九年六月

(2)

(3)

致謝

首先感謝指導教授邱奕鵬博士這兩年多來的教育，無論是在學術上或是研究態度上，都使我獲益良多，並給予我許多機會接觸學習新的事物，同時也要感謝口試委員林晃巖博士與王子建博士在口試當天針對論文內容給予許多相當有益的建議和指導。

另外，實驗室愉快的氣氛也是我碩士生涯中不可缺少的重要動力。感謝文嵐學姐、承翰學長與俊宏學長在研究上的協助，戰友們凡城、逸民、俊霖和建鴻學長帶給我的許多歡笑，以及早一步畢業的信毅學長、新生代禹秀學妹和牧珉學弟為實驗室增添的更多熱鬧和趣味。這一切，都將讓我難以忘懷。

感謝家中可愛的姊姊和妹妹，讓我的窩總是如此溫暖，話題也總是無窮無盡，

也謝謝同時在生理所打拼的曉君，能與我一路上相互砥礪。最重要的，是要謝謝我最愛的爸媽，除了生活上的照顧與關懷之外，也是我心中無比堅固的支柱和靠山，有他們陪伴是我生命中最大的驕傲。

(4)

摘要

本論文提出步階式數學形式之有限差分模態法來對光柵結構進行模擬並將其與傅立葉模態法(亦稱嚴格耦合波分析法)進行比較。我們的數值結果證實，對於橫電極化的所有情況以及橫磁極化在高導電性與無損金屬材料中，有限差分模態法比起嚴格耦合波分析法會有更好的收斂性和準確性。

對於有限差分模態法，我們考慮任意高階之邊界條件並將其與泰勒展開式結合。在不使所取格點數增加的情況下，我們亦將廣義道格拉斯(Douglas)方法套入使用來加速誤差收斂階數。使用前述技巧，可以對結構的每一層架構出稀疏矩陣並計算出存在於該層之模態所對應的場值分佈以及傳播常數。另外，我們也使用穆哈拉姆(Moharam)所提出之改良穿透矩陣方法來穩定多層光柵或甚至單層光柵層間的矩陣運算。

為了評估此種數值方法的可用性，我們將討論一些光柵的繞射特性，如入射角變化、厚度變化、佔空比變化所造成的影響，以及準確性、收斂性等等。另外，

我們也使用結合週期性邊界和吸收邊界的二維有限差分法來與前述方法進行比較。

關鍵詞：

傅立葉模態法、嚴格耦合波分析法、頻域有限差分法、任意階數邊界條件、廣義道格拉斯方法、改良穿透矩陣方法、完美匹配層。

(5)

Abstract

In this thesis, the finite-difference modal method (FDMM) with step- index formulation for simulating grating structures is proposed and compared with rigorous coupled-wave analysis (RCWA), also called Fourier modal method (FMM). It is verified that FDMM has better convergence and accuracy than RCWA for TE polarization in almost all cases and TM polarization for high conductive and lossless metallic materials.

In the FDMM, the relation of interface conditions to arbitrary high orders is considered and combines with Taylor series expansion. The generalized Douglas (GD) scheme is also adopted to enhance the convergence order without considering more sampled points. With the techniques mentioned above, the sparse matrix of eigenvalue problem could be constructed to solve the fields and the propagation constants of modes inside each layer. In addition, the enhanced transmittance matrix approach proposed by Moharam et al. for RCWA is used to make matrix manipulation stable for multi-layer or even single layer gratings.

The diffraction properties of gratings, such as accuracy, convergence, dependence of diffraction efficiencies on incident angle, thickness, duty cycle, etc, will be discussed for numerical assessment of FDMM. Moreover, two- dimensional finite-difference methods combined with periodic boundary conditions and absorbing boundary conditions will be executed for comparison.

(6)

Keywords:

Fourier modal method, rigorous coupled-wave analysis, finite-difference frequency-domain method, arbitrary-order interface conditions, generalized Douglas scheme, enhanced transmittance matrix approach, perfectly matched layer.

(7)

List of Figures

2.1 Grating with one-dimensional periodicity. . . 8

2.2 Rectangular-groove grating. . . 10

2.3 Example of multi-layer approximation. . . 16

2.4 One-dimensional discretization. . . 22

2.5 Sketch of a discontinuity between sampled points. . . 23

2.6 Periodic boundaries. . . 27

2.7 Configuration of two-dimensional finite difference. . . 31

2.8 Two-dimensional discretization. . . 32

2.9 Method of averaging permittivity. . . 34

2.10 Continuity of fields on interfaces. . . 35

2.11 A real region surrounded by PML to approximate the infinite region. . . 39

2.12 Configuration of multi-stretched-coordinate in PML. . . 43

2.13 Total-field and scattered-field regions connected by a plane wave source. . . 45

(10)

3.1 Fig.8 in the paper of Sheng et al. [30]. The zeroth-order reflection with respect to the incident angle. . . 51 3.2 Incident angle dependence of the zeroth-order reflection for

lossy metal with TM polarization (compared with Fig. 3.1). . 51 3.3 Field diagram of Fig. 3.2 at incident angle θ = 0^◦ (Left:

RCWA, Right: FDMM). . . 52 3.4 Field diagram of Fig. 3.2 at incident angle θ = 15^◦ (Left:

RCWA, Right: FDMM). . . 52 3.5 Fig.9 in the paper of Pai and Awada [25]. Transmittance with

respect to the incident angle. . . 53 3.6 Incident angle dependence of the zeroth- and first-order trans-

mittance for the dielectric triangular grating under TE incidence (comparing with Fig. 3.5). . . 53 3.7 Fig.6 of paper of Sheng et al. [30]. The first-order reflection of

Al and Ag gratings with respect to the thickness. . . 55 3.8 Thickness dependence of first-order reflection for Al and Ag

gratings under TM incidence (compared with Fig. 3.7). . . 55 3.9 Fig.6 of paper of Sheng et al. [30]. The absorption of Al and

Ag gratings with respect to the thickness. . . 56 3.10 Thickness dependence of absorption of Al and Ag gratings

under TM incidence (compared with Fig. 3.9). . . 56

(11)

3.11 Fig.1 of paper of Lyndin et al. [19]. The minus-first-order reflection with respect to groove width. . . 57 3.12 Groove width dependence of minus-first-order reflection (com-

paring with Fig. 3.11). . . 58 3.13 Accuracy and convergence of β of fundamental mode. A di-

electric grating with εr,g1 = 3² under TE incidence. . . 63 3.14 Convergence of the zeroth-order reflection. A dielectric grating

with ε_r,g1 = 3² under TE incidence. . . 63 3.15 Accuracy and convergence of β of fundamental mode. A loss-

less metallic grating with εr,g1 = (−10i)² under TE incidence. . 64 3.16 Convergence of the zeroth-order reflection. A lossless metallic

grating with ε_r,g1= (−10i)² under TE incidence. . . 64 3.17 Accuracy and convergence of β of fundamental mode. A lossy

grating with ε_r,g1= (3.18 − 4.41i)² under TE incidence. . . 65 3.18 Convergence of the zeroth-order reflection. A lossy grating

with ε_r,g1 = (3.18 − 4.41i)² under TE incidence. . . 65 3.19 Accuracy and convergence of β of fundamental mode. A high

conductive grating metallic with ε_r,g1 = (1 − 40i)² under TE incidence and nonuniform discretization. . . 66

(12)

3.20 Convergence of the zeroth-order reflection. A high conductive metallic grating with ε_r,g1 = (1 − 40i)² under TE incidence and nonuniform discretization. . . 66 3.21 Configuration of triangular gratings. . . 67 3.22 Revised discretization for multi-layer approximation. . . 68 3.23 Convergence of the zeroth-order reflection. A dielectric trian-

gular grating with ε_g1= 3² under TE incidence. . . 70 3.24 Convergence of the zeroth-order reflection with proper dis-

cretization. A dielectric triangular grating with ε_g1 = 3²under TE incidence. . . 70 3.25 Convergence of the zeroth-order reflection with proper dis-

cretization. A lossless metal triangular grating with ε_g1 = (−10i)² under TE incidence. . . 71 3.26 Convergence of the zeroth-order reflection with proper dis-

cretization. A lossy triangular grating with ε_g1 = (3.18 − 4.41i)² under TE incidence. . . 71 3.27 Accuracy and convergence of β of fundamental mode. A di-

electric grating with ε_r,g1 = 3² under TM incidence. Grids are distributed by ratio of width. . . 74

(13)

3.28 Convergence of the zeroth-order reflection. A dielectric grating with ε_r,g1 = 3² under TM incidence. Grids are distributed by ratio of width. . . 74 3.29 Accuracy and convergence of β of fundamental mode. A loss-

less metallic grating with ε_r,g1 = (−10i)² under TM incidence.

Grids are distributed by ratio of width. . . 75 3.30 Convergence of the zeroth-order reflection. A lossless metallic

grating with ε_r,g1 = (−10i)² under TM incidence. Grids are distributed by ratio of width. . . 75 3.31 Accuracy and convergence of β of fundamental mode. A lossy

grating with ε_r,g1= (3.18 − 4.41i)² under TM incidence. Grids are distributed by ratio of width. . . 76 3.32 Convergence of the zeroth-order reflection. A lossy grating

with ε_r,g1 = (3.18 − 4.41i)² under TM incidence. Grids are distributed by ratio of width. . . 76 3.33 Accuracy and convergence of β of fundamental mode. A lossy

grating with ε_g1 = (0.22 − 6.71i)² under TM incidence. Grids are distributed by ratio of width. . . 77 3.34 Convergence of the zeroth-order reflection. A lossy grating

with ε_g1 = (0.22 − 6.71i)² under TM incidence. Grids are distributed by ratio of width. . . 77

(14)

3.35 Accuracy and convergence of β of fundamental mode. A highly conductive grating with ε_g1= (1 − 40i)² under TM incidence.

Grids are distributed by ratio of width. . . 78 3.36 Convergence of the zeroth-order reflection. A highly conduc-

tive grating with ε_g1 = (1 − 40i)² under TM incidence. Grids are distributed by ratio of width. . . 78 3.37 Continuity of E_x fields for TM polarization. . . 79 3.38 Convergence of the zeroth-order transmittance. A dielectric

triangular grating with ε_g1= 3² under TM incidence. . . 84 3.39 Convergence of the zeroth-order transmittance with proper

discretization. A dielectric triangular grating with ε_g1 = 3² under TM incidence. . . 84 3.40 Convergence of the zeroth-order reflection with proper dis-

cretization. A lossless metal triangular grating with ε_g1 = (−10i)² under TM incidence. . . 85 3.41 Convergence of the zeroth-order reflection with proper dis-

cretization. A lossy triangular grating with ε_g1 = (3.18 − 4.41i)² under TM incidence. . . 85 3.42 Duty cycle variation of the minus-first-order reflection. A di-

electric grating with ε_g1= 3²under TM incidence. 2DFD with averaging permittivity is used. . . 86

(15)

3.43 Duty cycle variation of the minus-first-order reflection. A lossless metallic grating with ε_g1 = (−10i)² under TM incidence.

2DFD with averaging permittivity is used. . . 87 3.44 Duty cycle variation of the minus-first-order reflection. A loss-

less metallic grating with ε_g1 = (−10i)² under TM incidence.

2DFD with considering boundary condition is used. . . 88 3.45 Field diagram of Fig. 3.12 at groove width= 302nm (Left:

RCWA, Right: FDMM) . . . 89 3.46 Field diagram of Fig. 3.12 at groove width= 302nm (Left:

Graded-index 2DFD, Right: Step-index 2DFD) . . . 89 3.47 Field diagram of Fig. 3.12 at groove width= 250nm (Left:

RCWA, Right: FDMM) . . . 90 3.48 Field diagram of Fig. 3.12 at groove width= 250nm (Left:

Graded-index 2DFD, Right: Step-index 2DFD) . . . 90

(16)

Chapter 1 Introduction

1.1 Overview

Diffraction of optical electromagnetic radiation by periodic structures is important for many engineering applications. Grating diffraction is central in the fields of integrated optics, holography, optical data processing, spectral analysis, etc. There are numerous numerical methods with variety of possible assumptions to analyze diffraction properties of gratings, and most of them are not only tools for solving mathematical equations but also imply some physical insights of the problems. Thanks to advancement of computer tech- nology, numerical methods and computer-aided design (CAD) become more important and convenient for finding optimized parameters for high cost or complicated experiments.

(17)

1.2 Literature Survey

The Fourier modal method (FMM) which is a kind of differential methods [1]

and often referred to as rigorous coupled-wave analysis (RCWA) was first formulated for planar gratings by Moharam and Gaylord [2] [3] [4], and then extended to surface-relief gratings [7] [6] and crossed-grating structures. It provides the exact solutions whose accuracy depends solely on the numbers of terms retained in the space-harmonic expansions of the fields.

Before this method was proposed, the most common differential methods were the coupled-wave approach [8] and the modal approach [9]. The coupled- wave approach, which expanding the solution into plane-wave components, had been known to offer a relatively simple formulation and superior physical insight into wave-diffraction phenomena. In this approach, several assumptions were made in order to obtain solutions such as neglecting the second order derivatives of the field amplitudes and retaining only one diffracted wave. The modal approach, which through eigenmode expansion, was a rigorous exact analysis but complicated mathematically. However, Magnusson and Gaylord [10] had shown that these two approaches are equivalent and the coupled-wave approach could become rigorous by including all diffracted waves in the formulation together with retaining the second derivatives of the electromagnetic fields. Therefore, RCWA was created by considering the coupled-wave approach without the assumptions as mentioned above.

(18)

As matching the boundary conditions between every layers and finding diffraction efficiencies, some instabilities would be introduced because taking the inverse of some ill-condition matrices. Several methods such as R-matrix, S-matrix, and other approaches which were compared systematically by Li [11] have been proposed to produce stable implementation for this problem, and the method proposed by Moharam et al. [12] in 1995 will be mentioned in section 2.1.3.

Another problem of RCWA discovered by Li and Haggans [13] is poor convergence as dealing with metallic gratings under TM incidence. The convergence was dramatically improved by the reformulation by Lalanne and Morris [14] and Granet and Guizal [15], and this improvement was explained by Li [16] mathematically. He proved the origin of poor convergence comes from the mistakes of using Fourier factorization and proposed the inverse rule of factorization to uniformly satisfy the boundary conditions in the grating region.

However, even using the modified factorization mentioned above, Popov et al. [17] discovered a numerical instability problem in this differential theory as applied to metallic gratings with very high conductivity under TM incidence.

They attributed these instabilities to the imperfect condition of matrices gen- erated by Fourier coefficients of permittivity distribution to be inverted. This statement was later questioned by Watanabe [18]. Some heuristical solutions

(19)

proposed by them were introducing artificial metal losses in order to damp the instabilities or applying two-step truncation [17]. Such strategy obvi- ously treats a neighboring but different electromagnetic problem. Lyndin et al. [19] established the link between the instabilities and the spurious modes corresponding to instable high order eigenvalues. They proposed a procedure of identification and filtration of these spurious modes, but such tracking of those artifactual modes are complicated. Furthermore, Guizal et al. [20] applied the reformulation of FMM with adaptive spatial resolution proposed by Granet [21] to approach the problem of highly conductive gratings under TM incidence and get even more stable solution.

Finite difference (FD) methods for solving partial difference are also used in electromagnetism for solving Maxwell’s equations. This method is not widely used in grating theory but often used to study the diffraction by aperiodic objects of finite dimension because of their suitability for incor- porating absorbing boundary conditions to limit the computational domain.

Generally, standard FD methods require a two-dimensional (2D) mesh for the discretization of a one-dimensional (1D) grating. After adding periodic boundary conditions, absorbing boundaries and sources, the fields of scatter- ing problems could be solved and used to find diffraction efficiencies.

However, Lalanne and Hugonin [22] proposed a very simple method for the analysis of 1D lamellar gratings. This method is solving eigenmodes in-

(20)

side the gratings region such as FMM (RCWA) but using finite-difference approach. And it is similar to numerical techniques that are based on finite- difference modal approaches and used in waveguide theories. They used a first-order method with averaging permittivity and found that their method is much inferior to the RCWA for dielectric gratings. In contrast, their method compared favorably with the RCWA for metallic gratings operating in the infrared regions of the spectrum, especially for TM polarization. They also proposed three crucial methods to accelerate the convergence: (1) proper interpolation for averaging permittivity, (2) mesh points on discontinuities and (3) non-uniform sampling near the discontinuities. Their numerical results indicated that FD approach offers rather good performance for highly conducting gratings and TM polarization.

1.3 Motivation

In the paper of Lalanne and Hugonin [22], they presented the finite-difference modal method and declared good performance for metallic gratings and TM polarization but worse convergence than RCWA for TE polarization. Inside the grating layer, what they used is an interpolation scheme that locally averages the permittivity, and they said the method of interpolation has a drastic impact on the convergence performance. In addition, they only considered first-order finite-difference and expected that faster convergence

(21)

rates can be achieved by using higher order method but the computational efficiency will be reduced because of increasing non-zero value in the eigen matrices.

Therefore, the finite-difference modal method (FDMM) with even higher order formulation and considering boundary conditions instead of interpolation is introduced and tested in this thesis. With the generalized Douglas scheme [23] [24], the convergence order can be increased without adding the mesh points, so the computational time and the computer memory does not be increased. For comparison, the two-dimensional finite-difference both for graded-index approximation and step-index approximation will be investi- gated.

1.4 Chapter Outline

There are three chapters following this introduction.

In chapter 2, the Fourier modal method (FMM), or rigorous coupled wave analysis (RCWA) is demonstrated first both for rectangular-groove gratings and surface-relief gratings. Surface-relief gratings, which also called arbitrary profile gratings, can use many layers of rectangular-groove gratings with different duty cycles to approximate. A stable approach will then be mentioned to solve the unstable problems as using above approximation.

Next, the finite-difference modal method (FDMM) will be proposed by us-

(22)

ing step-index formulation with or without generalized Douglas scheme inside each grating layer and stable approach used in RCWA. Last, the two- dimensional finite-difference (2DFD) with averaging permittivity and considering boundary condition will be applied to the same configuration which solved by RCWA and FDMM above. In addition, two different methods of perfectly matched layers and total fields/scattered fields (TF/SF) for adding sources will be mentioned.

In chapter 3, numerical results are given to assess the formulations men- tion in chapter 2. First, some papers’ results will be used to verify the correctness of FDMM and RCWA, and the problems of RCWA will be seen for lossless metallic or high conductive gratings under TM incidence. Second, accuracy and convergence of propagation constants and diffraction efficiencies of these two numerical methods will be defined and discussed for every kinds of materials, both of TE and TM polarization and both of rectangular- groove gratings and arbitrary profile gratings. Finally, 2DFD will be used to solve the same problem which suffers from serious instabilities with lossless metallic gratings under TM incidence.

Chapter 4 concludes this thesis.

(23)

Chapter 2 Formulation

The geometry of the problem we deal with is one dimensional periodic grating depicted in Fig. 2.1, which is separated into reflection region, grating region and transmission region. The normal to the interface between any adjacent two regions is along z direction, and the grating is periodic along x direction and infinite along y direction.

Figure 2.1: Grating with one-dimensional periodicity.

(24)

2.1 Fourier Modal Method (or RCWA)

The main idea of RCWA is to express both permittivity and electromagnetic fields by Fourier bases. In the grating region, the periodic relative permittivity (x) is expandable in a Fourier series of the form

(x) =X

h

_he^j^2π^Λ^hx, (2.1)

where Λ is grating period and h is the hth Fourier coefficient of the relative permittivity which can be found by

_h = 1 Λ

Z xi+Λ xi

(x)e^−j^2π^Λ^hxdx, (2.2)

where x_i is any initial position in the integration. Regarding the electromagnetic fields, it assumes that these could be expanded by Fourier bases along the periodic direction x and form a set of modes along the propagation direction z and are expressed as

ψ_m(x, z) ∼ X

i

C_mie^−jk^xi^x

!

× e^±jβ^m^z, (2.3)

where m is the index of mode, β_m is the propagation constant of mode m, k_xi = k₀sin(θ_incidence) + i^2π_Λ and C_mi is the contribution of mode m to the ith Fourier order.

After using above approximation to rewrite Maxwell’s equations, an eigenvalue problem is constructed and used to find every eigenmodes inside the

(25)

grating layer. Then we match boundary condition of tangential field components at each interface between two adjacent layers and calculate the diffraction efficiencies finally. More details of RCWA will be shown in the following subsections.

2.1.1 Planar Diffraction of Rectangular-Groove Grat- ings

For simplicity, only TM polarization (H_y,E_x,E_z), which means the magnetic field is perpendicular to the plane of incidence (which means x-z plane here), is going to be demonstrated, and the mathematical derivation of TE polarization is similar to that of TM. The structure of a rectangular-groove grating is depicted in Fig. 2.2. Assume the incidence wave is a plane wave given by

Figure 2.2: Rectangular-groove grating.

H_y^inc= e^−jk⁰^√^r,inc(sin θx+cos θz), (2.4)

(26)

where k₀ is the wave vector in free space, and θ is the incident angle. However, it also could be any other kinds of waves. And the solutions form of reflection and transmission are given by

H_y^R=X

i

Rie^−j(k^xi^x−k^R^zi^z) (2.5)

and

H_y^T =X

i

T_ie^−j[k^xi^x+k^T^zi^(z−t)], (2.6) where k_zi^R = pk²₀_r,inc− k²_xi, k_zi^T = pk₀²_r,tra− k²_xi and t is the thickness of grating. The electric fields in these regions can be obtained from Maxwell’s equation E = _jω¹

0r∇ × H.

In the grating region, we would like to solve modes propagated along direction z. Being analogous to (2.3), the fields may be expressed as

H_y^G=X

i

U_yi(z)e^−jk^xi^x, (2.7)

where U_yi(z) means the contribution of every modes to the ith space harmonic fields. After substituting it into Maxwell’s equations, we will get

E_x^G = αX

i

S_yi(z)e^−jk^xi^x, (2.8)

where α is normalization constant which will be decided later. For TM polarization, Maxwell’s equations E = _jω¹

0r∇ × H and H = _jωµ⁻¹

0∇ × E will be simplified to be

E_z^G= 1 jω₀_r(x)

∂H_y^G

∂x (2.9)

(27)

E_x^G = 1

jω₀_r(x) −∂H_y^G

∂z

!

(2.10)

−jωµ₀H_y^G = ∂E_x^G

∂z − ∂E_z^G

∂x . (2.11)

First considering (2.10). If we define 1

_r(x) =X

h

a_he^j^2π^Λ^hx (2.12)

and substitute (2.12), (2.7) and (2.8) into (2.10). we will obtain α · S_xi = −1

jω₀ X

p

a_ip∂U_yi

∂z , (2.13)

where aip ≡ a_p−i is Fourier coefficient of ⁻¹_r (x). Next, combining (2.11) with (2.9) and substituting (2.12), (2.7) and (2.8) into this combination, then we obtain

α · ∂S_xi

∂z = −jk_xi jω0

X

p

a_ip(−jk_xi)U_yp− iωµ₀X

i

U_yi. (2.14) Now we can see that if α ≡ _jω¹

0, the (2.13) and (2.14) will become simpler.

Finally, combining the (2.13) with (2.14) and eliminating S_xi, we obtain X

p

a_ip∂²U_yp

∂z² = k_xiX

p

a_ipk_xpU_yp− k₀²U_i, (2.15) or, in matrix form,

AU⁰⁰ = K_xAK_xU − k²₀U, (2.16) where A ≡ [aip] and Kx≡ diag[kxi]. Therefore, now we know the eigenvalue problem used to be solved is

U⁰⁰ = A⁻¹(K_xAK_xU − k₀²I)U. (2.17)

(28)

However, from the empirical advice of Lalanne and Morris [14] (by quasi- static limit description), G. Granet and B. Guizal [15] in 1996 and the mathematical foundation given by Li [16] that the matrix A in the parenthesis of (2.17) is better to be replaced by E⁻¹, where E ≡ [_ip], for improving convergence. This replacement will improve the convergence. Namely, we solve

U⁰⁰= A⁻¹(K_xE⁻¹K_xU − k₀²I)U ≡ −MU (2.18) finally instead of (2.17). If the field assumption is like (2.3), i.e., U_yi(z) ∼ P

mCmie^±jβ^m^z, we could know that U⁰⁰=

"

X

i

∂²U_yi(z)

∂z²

#

= −β²U, (2.19)

and be sure that the eigenvalue problem is MU = β²U, where β² is a diagonal matrix containing eigenvalues.

After solving the eigenvalue problem above, the space harmonics of tangential electric and magnetic fields are given by

U_yi =X

m

w_im{g_m⁺e^−jβ^m^z+ g_m⁻e^+jβ^m^(z−t^g⁾} S_yi = −X

p

X

m

a_ipw_pm(−jβ_m){g_m⁺e^−jβ^m^z− g_m⁻e^+jβ^m^(z−t^g⁾}, (2.20)

where w_im and β_m are the elements of the eigenvector matrix W and the square root of the eigenvalues of the matrix M. The quantities g_m⁺ and g_m⁻ are unknown constants to be determined by matching boundary conditions.

Physically, g⁺_m and g_m⁻ mean the contribution of each mode. The columns of

(29)

[w_im] means which mode, and the ith row of column m represents the ratio of tth space harmonic in this mode.

To solve reflection coefficient R_i, transmission coefficient T_i, g_m⁺ and g_m⁻, the boundary conditions of tangential fields are used. For TM polarization, the component E_x, which means ¹

r

∂Hy

∂z , and H_y are continuous at discontinuities. At the input boundary (z = 0)

δ_i0+ R_i =X

m

w_im{g⁺_m+ e^−jβ^m^t^gg_m⁻} (2.21) and

−jk0

√r,inccos θδi0

r,inc

+ jk^R_zi

r,inc

R_i =X

p

X

m

a_ipw_pm(−jβ_m){g_m⁺− e^−jβ^m^t^gg_m⁻}, (2.22) or in matrix form,

"

δ_i0

−jk√₀cos θ

r,inc δi0

# +

"

I

−Z_R

# R =

"

W WX

AWZ −AWZX

# "

g_m⁺ g_m⁻

#

, (2.23)

where Z_R ≡ diag[−jk_zi^R/_r,inc], Z ≡ diag[−jβ_m] and X ≡ diag[exp(−jβ_mt_g)].

At z = t_g,

X

m

w_im{e^−jβ^m^t^gg_m⁺+ g_m⁻} = T_i (2.24) and

X

p

X

m

aipwpm(−jβm){e^−jβ^m^t^gg_m⁺− g_m⁻} = −jk_zi^T

_r,tra Ti, (2.25) or in matrix form,

"

WX W

AWZX −AWZ

# "

g_m⁺ g_m⁻

#

=

"

I Z_T

#

T, (2.26)

(30)

where Z_T ≡ diag[−jk^T_zi/_r,tra].

The coefficients R_i and T_i could be found by using (2.23) and (2.26). The diffraction efficiencies could be calculated by the definition of the ratio of Poynting’s vectors along propagation direction z. For TM polarization here, the Poynting’s vector in the propagation direction z is

P_z = 1

2ReE_xH_y^∗ = 1 2Re

j

ω₀_r(∂H_y

∂z )H_y^∗

, (2.27)

and the diffraction efficiencies are defined as DEri = P_z^R

P_z^inc = |Ri|²Re k_zi^R k_z0

DEt_i = P_z^T

P_z^inc = |T_i|²Re _r,inc

r,tra

k_zi^T kz0

. (2.28)

For lossless gratings, the sum of the reflected and transmitted diffraction efficiencies given by (2.28) must be unity, which means conservation of energy.

In addition, if it needs to plot the field diagrams, solve coefficients g_m⁺ and g⁻_m after finding R_i and T_i. In the next subsection, we will generalize the structure to arbitrary shape.

2.1.2 Multi-layer Approximation for Arbitrary Shape Gratings

For gratings with arbitrary profiles, we divide the grating into a large number of sufficiently thin layers and approximate each layer by a rectangular-groove grating, which was use first proposed by Peng et al. and applied to RCWA by Moharam and Gaylord [6], as in Fig. 2.3.

(31)

Figure 2.3: Example of multi-layer approximation.

The electromagnetic fields in each grating layer are determined by RCWA for rectangular-groove gratings (or by other approaches, such as the modal approach and the finite difference approach). Then boundary conditions of the tangential fields are applied in sequence at every interfaces to get reflected and transmitted diffracted field amplitudes and diffraction efficiencies. Here we still use TM polarization to illustrate the formulation.

After solving the eigenvalue problems in each grating layer, the space harmonics of tangential fields of every layers are given by

U_l,yi =X

m

w_im^l {g_l,m⁺ e^−jβ^m^l ^(z−z^l−1⁾+ g_l,m⁻ e^+jβ^l^m^(z−z^l⁾} S_l,yi = −X

p

X

m

a^l_ipw_pm^l (−jβ_m^l ){g⁺_l,me^−jβ^m^l ^(z−z^l−1⁾− g_l,m⁻ e^+jβ^l^m^(z−z^l⁾}, (2.29)

where the index l means which grating layer. Then we match boundary con-

(32)

dition at every interfaces and express these in matrix form. At the interface between the input region and the first grating layer, i.e. at z = 0 in Fig. 2.3,

"

δi0

−jk₀cos θ

√r,inc δ_i0

# +

"

I

−Z_R

# R =

"

W₁ W₁X₁

A₁W₁Z₁ −A₁W₁Z₁X₁

# "

g_1,m⁺ g_1,m⁻

# , (2.30) where the suffixes of W, A, Z and X mean which grating layer and the definitions of these matrices are same as in preceding subsection, at the interface between (l − 1)th and lth layer (z = z_l−1)

"

Wl-1Xl-1 Wl-1

A_l-1W_l-1Z_l-1X_l-1 −A_l-1W_l-1Z_l-1

# "

g_l−1,m⁺ g_l−1,m⁻

#

=

"

W_l W_lX_l A_lW_lZ_l −A_lW_lZ_lX_l

# "

g_l,m⁺ g_l,m⁻

#

, (2.31)

and at the interface between the last grating layer and the output region

"

W_LX_L W_L

A_LW_LZ_LX_L −A_LW_LZ_L

# "

g⁺_L,m g⁻_L,m

#

=

"

I Z_T

#

T. (2.32)

Finally, (2.30)-(2.32) are rewritten as

"

δ_i0

−jk0cos θ

√r,inc δ_i0

# +

"

I

−Z_R

# R =

L

Y

l=1

"

W_l W_lX_l A_lW_lZ_l −A_lW_lZ_lX_l

#

×

"

WlXl Wl

A_lW_lZ_lX_l −A_lW_lZ_l

#⁻¹"

I Z_T

#

T. (2.33)

However, this approximation used in successive field matching may introduce numerical instability which is due to the presence of evanescent fields.

And these evanescent waves, which possess large imaginary part of β, will

(33)

make the matrix X to be ill-conditioned as doing matrix inverse. In the next subsection we will introduce stable approaches to eliminate the numerical instability.

2.1.3 Stable Approach for Multi-layer Approximation

There are several approaches which have been proposed to produce stable implementation. For example, Moharam and Gaylord [6] obtained numeri- cally stable RCWA calculation for TE polarization and dielectric gratings to a grating depth of as many as four wavelengths, by sequential Gaussian elim- ination scheme. Pai and Awada [25] use layer transmission matrices and interface reflection and transmission matrices to derive the solution for RCWA in terms of a multiple-reflection series which is stable for TE polarization and dielectric gratings to a grating depth of as many as four wavelengths. Li [26]

used the R-matrix propagation algorithm to propagate the field through the layers in the modal approach to obtain stable results for deep dielectric and metallic one-dimensional gratings in the conical mount.

Here we use the method which is proposed by Moharam, Pommet, and Grann [12] in 1995 and called enhanced transmittance matrix approach. See (2.33) again and notice the inverse matrix at left. There are some X terms in this big inverse matrix. Because the elements of X are e^±jβ^m^z, as some βm values have large imaginary part, these exponential terms will produce very large or very small elements in this big inverse matrix. And the inver-

(34)

sion of this almost singular matrix will produce erroneous results because of truncation errors. Therefore, we would like to remove X terms from the big inverse matrix.

Consider the last factor (l = L) in (2.33), which is

"

WL WLXL

A_LW_LZ_L −A_lW_lZ_lX_l

# "

WLXL WL

A_LW_LZ_LX_L −A_lW_LZ_L

#⁻¹"

fT

g_T

# T

=

"

W_L W_LX_L

ALWLZL −AlWlZlXl

#

×

"

X_L 0

0 I

#⁻¹"

W_L W_L

A_LW_LZ_L −A_lW_LZ_L

#⁻¹"

f_T g_T

#

T,(2.34)

where f_T = I and g_T = Z_T. The matrix to be inverted has been rewritten as the product of two matrices. The matrix on the right side in the product is well conditioned, however, the left side matrix is ill-conditioned. By defining

"

W_L W_L

A_LW_LZ_L −A_lW_LZ_L

#⁻¹"

f_L g_L

#

≡

"

u_L v_L

#

, (2.35)

(2.34) could be reduced to

"

W_L W_LX_L

ALWLZL −AlWlZlXl

# "

X⁻¹_L 0

0 I

# "

u_L vL

# T

=

"

W_L W_LX_L

# "

X⁻¹_L u_L v_L

#

T, (2.36)

and by changing variable T = u⁻¹_L X_LT_L, the (2.36) will become

"

W_L W_LX_L

# "

I v_Lu⁻¹_L X_L

# T_L ≡

"

f_L g_L

#

T_L. (2.37)

(35)

Repeat this process for all layers, we obtain

"

δ_i0

−jk√0cos θ

r,inc δ_i0

# +

"

I

−ZR

# R =

"

f₁ g1

#

T₁, (2.38)

where T = u⁻¹_L X_Lu⁻¹_L−1X_L−1. . . u⁻¹₂ X₂u⁻¹₁ X₁T₁. In this formulation, the instability could be avoided successfully because it never inverts the matrices X. The singular-value decomposition technique could be considered in inverting the matrix u to avoid numerical difficulties because of round-off errors when a large number of layers and a large number of harmonics are used.

In addition, if it needs to plot the field diagram, we have to find coefficients g⁺_l,m and g_l,m⁻ in every layers by finding R, T and Tl of each layer first and then substituting back. After using the skill above, the coefficients in each layer will be

"

g_L,m⁺ g_L,m⁻

#

=

"

I vLu⁻¹_L XL

# T_L,

"

g_L−1,m⁺ g_L−1,m⁻

#

=

"

I

vL−1u⁻¹_L−1XL−1

#

T_L−1, · · · (2.39) , where T_L = u⁻¹_L−1X_L−1· · · u⁻¹₁ X₁T₁, T_L−1 = u⁻¹_L−2X_L−2· · · u⁻¹₁ X₁T₁ and so on.

(36)

2.2 Finite-Difference Modal Method

In the finite difference modal method, we transfer the continuous functions and their differentiations into discrete points values. The assumption of fields is to express the distribution of spatial fields in grating period direction x and, just like RCWA method, form a set of modes in the propagation direction z which is

ψm(x, z) ∼ ψm(x) × e^±jβ^m^z, (2.40) where m is the index of mode, β_m is the propagation constant of mode m.

The equation to be solved is Helmholtz equation, ∇²ψ + k²₀_rψ = 0 (as- suming non-magnetic material). If the structure is infinite and uniform in y-axis and periodic in x-axis, and the direction of fields is in y-axis (i.e. TE or TM polarization), the equation becomes

∂²ψ_y(x, z)

∂x² +∂²ψ_y(x, z)

∂z² + k₀²_r(x)ψ_y(x, z) = 0. (2.41) Next, substituting (2.40) into (2.41), we obtain

∂²ψ_y(x)

∂x² + k₀²_r(x)ψ_y(x) = β²ψ_y(x). (2.42) After discretizing (2.42), constructing sparse matrix, and solving eigenvalue problem inside each layer, eigenmodes with x distribution inside each layer are found. Then boundary conditions are matched to solve the diffraction efficiencies and plot the field diagrams. In the next subsection, more details about the finite difference inside each layer are offered.

(37)

2.2.1 Eigenvalue Problems inside Each Layer

The basic idea of finite difference is to express the differentiation of a field point by itself and it’s adjacent points, as in Fig. 2.4. Being different from

Figure 2.4: One-dimensional discretization.

RCWA, the convergence of FD method could be improved by considering more adjacent points, i.e. extending to higher orders of Taylor series expansion. If we use (2N + 1) points to approximate the differentiation of one field point, ψ_i, we have to know the relation between the field at those points and ψ_i and its derivatives up to (2N )th order as





 ψ_i−N

... ψ_i

... ψ_i+N







=







u_−N,0 u_−N,1 . . . u_−N,2N ... ... . .. ... u_N,0 u_N,1 . . . u_N,2N











 ψ_i

... ψ_i^(j)

... ψ_i^2N







+ O(h^{2N +1}). (2.43)

And then inverting the matrix to find ψ_i⁰⁰with truncation error O(h^{2N +1}/h²) = O(h^{2N −1}). As the grids are uniformly positioned, the truncation error will become O(h^2N). Here we use central difference scheme, but it is not necessary.

We can use forward-difference or backward-difference as well.

(38)

Figure 2.5: Sketch of a discontinuity between sampled points.

However, there are some discontinuities in each grating layer. Unlike RCWA, which use Fourier expansion to illustrate the change of permittivity, here we match boundary conditions at these discontinuities to arbitrary higher order. As in Fig. 2.5, if we want to express the (i + 1)th points in terms of each order of the ith point, we expand each order of (i + 1)th field ψ_i+1 by each order of ψ_R, i.e.

Ψ_i+1≡





 ψ_i+1 ψ_i+1⁰ ψ_i+1⁰⁰ ψ_i+1⁽³⁾

... ψ_i+1^{(2N )}







=







1 q ^q_2!² ^q_3!³ · · · _{(2N )!}^q^2N 0 1 q ^q_2!² · · · _{(2N −1)!}^q^{2N −1} 0 0 1 q · · · _{(2N −2)!}^q^{2N −2} 0 0 0 1 · · · _{(2N −3)!}^q^{2N −3}

... ... ... ... . .. ... 0 0 0 0 · · · 1











 ψ_R ψ_R⁰ ψ_R⁰⁰ ψ_R⁽³⁾

... ψ^{(2N )}_R







≡ M_i+1:RΨ_R,

(2.44)

(39)

and each order of ψ_L by each order of (i)th field ψ_i,

Ψ_L ≡





 ψ_L ψ_L⁰ ψ_L⁰⁰ ψ_L⁽³⁾

... ψ_L^{(2N )}







=







1 p ^p_2!² ^p_3!³ · · · _{(2N )!}^p^2N 0 1 p ^p_2!² · · · _{(2N −1)!}^p^{2N −1} 0 0 1 p · · · _{(2N −2)!}^p^{2N −2} 0 0 0 1 · · · _{(2N −3)!}^p^{2N −3}

... ... ... ... . .. ... 0 0 0 0 · · · 1











 ψ_i ψ⁰_i ψ_i⁰⁰ ψ_i⁽³⁾

... ψ^{(2N )}_i







≡ M_L:iΨ_i.

(2.45) And the main problem is how to connect Ψ_L and Ψ_R, i.e., because

Ψ_i+1 = M_i+1:RM_R:LM_L:iΨ_i, (2.46)

we need to find the matrix MR:L. Next, we will derive the arbitrary higher order boundary condition for (2.42).

A. Arbitrary Higher Order Boundary Condition

[23] [24]

Generally, considering the discontinuous parameters of materials are not only permittivity but also permeability µ. The zeroth-order boundary condition is continuing of tangential fields

ψ_R= ψ_L, (2.47)

where ψ = E_y for TE and H_y for TM, and subscript R and L mean points infinitesimally close to the interface on the left and right, respectively. Next, for the first-order boundary condition, we use Maxwell’s equations H =

−1

jωµ0µr∇ × E and H_z,R= H_z,L for TE; E = _jω¹

0r∇ × H and E_z,R= E_z,L for

(40)

TM to derive

ψ_R⁰ = aψ⁰_L, (2.48)

where a = ^µ_µ^R

L for TE and a = ^R

L for TM. The second-order boundary condition is derived from Helmholtz equation, i.e. ψ_R⁰⁰ + ω²µ_R_Rψ_R = β²ψ_R and ψ⁰⁰_L+ ω²µ_L_Lψ_L = β²ψ_L. By using (2.47), we could obtain

ψ⁰⁰_R= ψ_L⁰⁰+ bψ_L, (2.49)

where b = ω²(µ_L_L− µ_R_R). For the third-order boundary condition, the Helmoholtz equation is differentiated on each side to get ψ_R³ + ω²µRRψ⁰_R= β²ψ⁰_R and ψ³_L+ ω²µ_L_Lψ⁰_L= β²ψ⁰_L. After substituting (2.48) into above two equations, we obtain

ψ_R⁽³⁾ = a(ψ_L⁽³⁾+ bψ_L⁰ ). (2.50) We proceed to find the 4th order boundary relation. In the derivation here, we could regard the term (∂_x²+ ω²µ) as an operator, which is equivalent to the term β², i.e.,