特殊光波導的分析與設計(IV)

行政院國家科學委員會專題研究計畫 成果報告

特殊光波導的分析與設計(IV)

計畫類別： 個別型計畫

計畫編號： NSC93-2215-E-002-042-

執行期間： 93 年 08 月 01 日至 94 年 10 月 31 日

執行單位： 國立臺灣大學光電工程學研究所

計畫主持人： 張宏鈞

計畫參與人員： 許森明(台大光電所) 周允賢(台大光電所) 黃介銘(台大光

電所) 江瑞倫(台大電信所) 陳明昀(台大光電所)

報告類型： 精簡報告

處理方式： 本計畫可公開查詢

中 華 民 國 95 年 4 月 18 日

1

行政院國家科學委員會專題研究計畫成果報告

特殊光波導的分析與設計(IV)

Analysis and Design of Special-Type Optical Waveguides (IV)

計畫編號：NSC 93-2215-E-002-042

執行期限：93 年 8 月 1 日至 94 年 10 月 31 日

主持人：張宏鈞 台灣大學電機系、光電所暨電信所教授

計畫參與人員：許森明(台大光電所) 周允賢(台大光電所)

黃介銘(台大光電所) 江瑞倫(台大電信所)

陳明昀(台大光電所)

一、中文摘要

本計畫將已建立的採用線性複合邊界

結點三角形元素的全向量有限元素法推展

為採用曲線性三角形元素的程式，並維持

完全匹配層吸收邊界，適合應用於如各式

光纖之具曲形材質介面的結構的更精準分

析。所建立的數值模型得以分析各式光子

晶體孔洞光纖的向量模態以及洩漏特性，

本計畫研究了不同空氣孔洞圈數的孔洞光

纖以及不同圈數的光子能隙光孅。本計畫

並發展具完全匹配層的二維有限元素波束

傳播法，得以模擬包括抗諧振反射光波導

之板形波導的各項特性。

關鍵詞：光波導、光子晶體光纖、抗諧振

反射光波導、有限元素法、波束

傳播法、光波導元件

Abstract

In this research the full-vectorial finite

element method (FEM) based on linear

hybrid edge/nodal elements with triangular

shape and having perfectly matched layer

(PML) absorbing boundaries will be

generalized by employing curvilinear

elements for more accurately analyzing

waveguides with curved boundaries, such as

optical fibers and photonic crystal fibers.

The established numerical model can analyze

vector modes and leaky characteristics of

various photonic crystal holey fibers. We

have studied holey fibers and photonic band

gap fibers of different numbers of air-hole

rings. A two-dimensional beam propagation

method model based on the finite element

method and with perfectly matched layer

absorbing boundaries has also been

developed for investigating various

properties of slab waveguides including

ARROWs.

Keywords: Optical waveguides, photonic

crystal fibers, antiresonant

reflecting optical waveguides

(ARROWs), finite element

method, beam propagation

method, optical waveguide

devices

二、計畫緣由與目的

本計畫以有限元素法為基礎，發展高

準度的波導模態分析數值模型，以輔助研

發及設計幾種特殊的光導波結構，主要的

探討對象為各式光子晶體光纖與孔洞光

纖。先前我們已建立了採用線性複合邊界

結點三角形元素的全向量有限元素法的數

值模型[1]，在本研究中我們進一步推展為

採用曲線性三角形元素的程式，並維持完

全匹配層吸收邊界，適合應用於如各式光

纖之具曲形材質介面的結構的更精準分

析。我們以此模型分析各式光子晶體孔洞

光纖的向量模態以及洩漏特性，並成功計

算不同空氣孔洞圈數的孔洞光纖以及不同

圈數的光子能隙光孅。除了特徵值問題求

解模型外，本計畫並發展具完全匹配層的

二維有限元素波束傳播法，得以另一方式

計算包括抗諧振反射光波導之板形波導的

模態特性。本報告根據曲線性三角形元素

2

有限元素分析法的建立、發展二維的具

PML 的有限元素波束傳播法、以及光子晶

體光纖的研究等三部分分別說明研究成

果。

三、結果與討論

　　

1. 曲線性三角形元素有限元素分析法的建

立

我們已建立有採用線性複合邊界結點

三角形元素的全向量有限元素法，本研究

中我們已推展為採用曲線性三角形元素的

程式，與[2]不同的是我們維持 PML 吸收邊

界。應用於如光纖之具曲形材質介面的結

構時，在波導截面先做線性複合邊界結點

三角形元素分割，然後將曲形材質介面鄰

近的元素變形為如圖一所示之曲線性三角

形元素，可增大程式的效率與精準度。由

於保有 PML 邊界，由推得的特徵值問題將

得以求得可能的複數傳播常數。圖二為強

導波光纖的例子，圖二(a)為結構截面圖，

外圍加了 PML，圖二(b)為四分之一區域的

有限元素分割分布圖，因結構對稱關係，

計算區域只須考慮此四分之一區域即可。

圖二(b)顯示線性有限元素分布，但在纖芯

與空氣纖殼接面的元素，我們會將三角形

元素的近接面一邊推成曲線形以更匹配曲

形接面。表一列出以全向量有限元素法模

態分析之基模的有效折射率 n

eff

= β/k

0

與正

規化傳播常數[(β/k

0

)

2

-n

22

]

2

/(n

12

-n

22

)隨有限

元素未知數數目的收斂情形，其中β為傳播

常數，k

0

為真空中的波數。表一顯示，與

正確值比較，本有限元素法分析之正規化

傳播常數計算有六位準確度。

2. 發展二維的具 PML 的有限元素波束傳

播法

與有限差分波束傳播法相同，波動方

程式經由近軸假設並引進緩慢變化脈波波

形 近 似 (slowly varying envelope

approxi-mation)，可導出沿傳播方向的一階微分方

程式。垂直於傳播方向的二階微分項，在

有限差分波束傳播法中係以數值為差分處

理。在有限元素波束傳播法中，仍以有限

差分法求解沿傳播方向的一階微分方程

式，但垂直於傳播方向的問題則以有限元

素法處理，其元素分割與上述模態分析模

型相同。我們發展二維的具 PML 的有限元

素波束傳播法，並建立虛軸波束傳播法

(imaginary-distance BPM, or ID-BPM) [3]，

得以模擬二維板形洩漏型波導，如二維抗

諧 振 光 波 導 (ARROW) 的 各 項 特 性 。

ARROW 的傳播常數為複數，代表光波沿

波導具有功率損耗，近年來受到許多注

意，它的優點包括有效的單模態傳播、低

損耗、較大的模場分佈適於與光纖連接、

容易製造、具有用的極化及波長選擇特性

等[4], [5]。首先檢驗一板形洩漏型波導之

計算，圖三為其折射率分布圖，圖四為以

有限元素模態解法所得到的基模場型圖，

有效折射率的實部為 1.4462244197，虛部

為 0.0004721871。圖五為以二維的具 PML

的有限元素波束傳播法模擬此基模的演進

情形，表二列出有效折射率的實部與虛部

隨傳播步驟的收斂情形，顯示與模態解法

之結果非常吻合。圖六為一二維 ARROW

的折射率分布圖，其前四個 TM 模態的場

型示於圖七。

3. 光子晶體光纖的研究

我們應用所建立的全向量模態分析模

型計算數種光子晶體孔洞光纖的向量模態

傳播特性以及洩漏特性，並成功計算不同

空氣孔洞圈數的孔洞光纖以及不同圈數的

光子能隙光孅。詳細內容請參附錄一[6]。

四、計畫成果自評

本計畫研究內容與原計畫相符，預期

目標大致達成。本計畫將原已建立的採用

線性複合邊界結點三角形元素的全向量有

限元素法推展為採用曲線性三角形元素的

程式，並維持完全匹配層吸收邊界，因而

對於如各式光纖之具曲形材質介面的結構

得以獲得更精準分析。此一數值模型可應

用於分析各式光子晶體孔洞光纖的向量模

態以及洩漏特性，本計畫展示了不同空氣

孔洞圈數的孔洞光纖以及不同圈數的光子

能 隙 光 纖 的 分 析 成 績 ， 已 發 表 於 Proc.

SPIE。本計畫並發展具完全匹配層的二維

有限元素波束傳播法，包括虛軸波束傳播

法，得以非常準確地計算包括抗諧振反射

3

光波導之板形波導的各項特性，如洩漏特

性。

五、參考文獻

[1]

H. J. Chen and H. C. Chang, “Modal

analysis of two-dimensional antiresonant

reflecting optical waveguides,” in

Proceedings of Optics and Photonics

Taiwan ’03 (OPT ’03), Vol. III, pp. 212–

214, Taipei, Taiwan, R.O.C., December

25–26, 2003.

[2]

M. Koshiba and Y. Tsuji, “Curvilinear

hybrid edge/nodal elements with

triangular shape for guided-wave

problems,” J. Lightwave Technol., vol.

18, pp. 737-743, 2000.

[3]

M. M. Spühler, D. Wiesmann, P. Freuler,

and M. Diergardt, “Direct computation

of higher-order propagation modes using

the imaginary-distance beam

propagation method,” Opt. Quamtum

Electron., vol. 31, pp. 751–761, 1999.

[4]

M. A. Duguay, Y. Kokubun, T. L. Koch,

and L. Pfeiffer, “Antiresonant reflecting

optical waveguides in SiO

2

-Si multilayer

structures,” Appl. Phys. Lett., vol. 49, pp.

13-15, 1986.

[5]

T. Baba and Y. Kokubun, “Dispersion

and radiation loss characteristics of

antiresonant reflecting optical

waveguides--Numerical results and

analytical expressions,” IEEE J.

Quantum Electron., vol. 28, pp.

1689-1700, 1992.

[6]

S. M. Hsu, H. J. Chen, and H. C. Chang,

“Finite Element Analysis of

Full-Vectorial Modal and Leakage properties

of Microstructured and Photonic Crystal

Fibers,” Proceedings of SPIE, Vol. 5623,

pp. 316–324, 2005.

六、圖表

表一

Effective index and normalized propagation constant calculation using FEM solver

for a strongly guiding fiber

0.82187881

1.43684306

exact

0.821878526690

121

1.436842933215

81

15101

0.821877984

217860

1.436842688

71320

9457

0.821875888240

527

1.4368417440157

9

4509

0.821859840182

287

1.436834510825

83

2897

0.8218540125907

55

1.4368318842014

4

2753

b

n

eff

unknowns

4

表二

Numerical convergence of real and imaginary parts of the effective index for the

ARROW of Fig. 3 calculated using the ID-FE-BPM

圖一 曲線性複合邊界結點三角形元素[2]。

50

49

48

7

6

5

4

3

2

1

step

0.00047218711051399

1.44622441977636

0.000472187110513987

1.44622441977636

0.00047218711051398

1.44622441977636

0.000473912073562678

1.44622457847957

0.000476342878347758

1.44622409797834

0.000477946578256607

1.4462266187241

0.000486920578214274

1.44621716690528

0.000493907731252267

1.4462488724262

0.000471868364156626

1.44613245563229

0.000665231175849808

1.44627486187387

Im[n

eff

]

Re[n

eff

]

5

圖二 強導波光纖的例子。(a)結構截面圖，外圍加了 PML。(b)四分之一區域的有限元素

分割分布圖。

圖三 板形洩漏型波導之折射率分布圖。

圖四 以有限元素模態解法所得到的圖三波導基模場型圖

-4 -3 -2 -1 0 1 2 3 4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 y A m p lit u d e

6

圖五 以二維的具 PML 的有限元素波束傳播法模擬圖三波導基模的演進情形。

圖六 二維 ARROW 的折射率分布圖。

TM

1

TM

2

TM

3

TM

4

圖七 圖六之二維抗諧振光波導前四個 TM 模態的場型圖。

-4 -3 -2 -1 0 1 2 3 4 0 10 20 30 40 50 60 70 80 y ste p

Finite element analysis of full-vectorial modal and leakage properties

of microstructured and photonic crystal fibers

Sen-Ming Hsu, Hung-Jen Chen, and Hung-chun Chang

Department of Electrical Engineering, Graduate Institute of Electro-Optical Engineering, and

Graduate Institute of Communication Engineering

National Taiwan University, Taipei, Taiwan 106-17, Republic of China

E-mail:

hcchang@cc.ee.ntu.edu.tw

ABSTRACT

Full-vectorial finite element method based eigenmode solver and imaginary-distance beam propagation method are developed for analyzing modal leakage characteristics of microstructured fibers and photonic band gap fibers. In both schemes the curvilinear hybrid edge/nodal elements based on linear tangential and quadratic normal vector basis functions are adopted to accomplish the computational window divisions and perfectly matched layer (PML) is incorporated as the boundary condition to absorb waves out of the computational window. The two schemes give consistent results for the numerical examples considered and the validity and usefulness of this work are demonstrated. Comparison with published results is shown.

Keywords: Finite element method, beam propagation method, perfectly matched layers, imaginary distance propagation method, photonic crystal fibers, microstructured fibers, holey fibers, photonic band gap fibers, optical waveguides, electromagnetic theory

1. INTRODUCTION

Microstructured fibers (MFs) and phtonic crystal fibers (PCFs) have been rapidly attracting intensive research interest in fiber optics.1,2 The MFs are featured by several air holes running parallel to the fiber axis in their structures such that the holes serve as low-refractive-index regions (lower than the index of silica or other polymer materials) to make an equivalent lower-index cladding compared to the central guiding region. The number of air holes and the arrangement of the holes determine the propagation characteristics of such fibers. The PCFs may have their cladding composed of a two-dimensional periodic distribution of air holes. If the light guiding in the central “defect” core region is caused by the photonic band gap (PBG) effect of the cladding, they are named as PBG fibers.3 It is now well known that the propagating modes of these novel fibers possess confinement or leaky losses, or their modal propagation constants are complex.4 Among various proposed analysis methods for such fibers, the finite element method (FEM) based ones have proved to be very useful and of high numerical accuracy.5,6 In this paper, we report our development of the FEM codes for the analysis of the MFs and the PBG fibers.

Propagating modes of optical waveguides can be calculated by either the eigenmode solver or the beam propagation method (BPM).7 The eigenmode solver solves the eigenvalue problems associated with the wave equations, while the BPM calculates the modes by a propagating scheme. In Ref. 5, a full-vectorial imaginary-distance BPM under the finite element formulation, called the FE-ID-BPM, was developed for calculating the real and imaginary parts of the propagation constants. In Ref. 6, an eigenvalue problem was formulated through the FEM and the complex eigenvalues were solved by a particular technique.8 We have established FEM-based eigenmode solver and BPM codes. In the derivation of the FEM-based eigenmode solver, we restrict the permittivity and permeability tensors to the diagonal form to obtain a simple eigenvalue equation and solve it. On the other hand, in the formulation of the FEM-based BPM, we assume the material tensors to be non-diagonal for more general cases and utilize the concept of the imaginary-distance propagation method7 to form the FE-ID-BPM. Using the two developed methods, we analyze MFs and PBG fibers and calculate their complex modal propagation constants. The results obtained using the two methods are compared. We also compare our results with those of Refs. 5 and 6, and some slight difference from Ref. 6 for the PBG fibers is observed.

Passive Components and Fiber-based Devices, edited by Yan Sun, Shuisheng Jian, Sang Bae Lee, Katsunari Okamoto, Proc. of SPIE Vol. 5623 (SPIE, Bellingham, WA, 2005) · 0277-786X/05/$15 · doi: 10.1117/12.576116 316

In both the ID-BPM and the eigenmode solver, the fiber waveguide cross-section is divided into curvilinear hybrid edge/nodal elements with triangular shape. Such elements are particularly useful for imposing the continuity of the tangential fields. Linear tangential and quadratic normal (LT/QN) basis functions are employed. The curved boundaries are modeled by the curvilinear elements for achieving better accuracy and efficiency. The adaptive mesh generation tool with Delaunay triangulation, named EasyMesh developed by B. Niceno,9 is employed in this work. For obtaining leaky losses, perfectly matched layers (PMLs) for anisotropic media10 are incorporated as the absorbing boundary of the computational domain. In the eigenmode solver, we employ the shift inverse power method (SIPM) for searching for the eigen solutions. However, for the current problem, since the desired mode is not the fundamental mode with the largest mode index, the SIPM would converge to a wrong mode index if we used the refractive index of the fiber material (say, silica) as the initial guessing index for iteration. We have made a detailed observation on the convergence characteristics of the mode-searching processes based on the SIPM and modified the SIPM so that we can solve the desired eigenmode of the leaky waveguide problem more efficiently. The desired mode would have mode field distributed around the center of the fiber. After several modified steps the correct mode index could be reached. The numerical formulations are described in Section 2 and numerical results are given in Section 3, followed by the conclusion in Section 4.

2. FORMULATION

2.1. Basics of the finite element eigenmode solver

Referring to Fig. 1 and the related coordinate system, we consider an anisotropic optical waveguide with its transverse material properties described by the relative permittivity and the permeability tensors

[ ]

0 0 00 0 0 rx r ry rz ε ε ε ε     =_{ }     ,

[ ]

0 0 0 0 . 0 0 rx r ry rz µ µ µ µ     =_{ }     (1)

Using PMLs for anisotropic media,10 we have the vectorial wave equation derived from Maxwell’s equations:

[ ]

(

)

2

[ ]

0 0

p k q

∇× ∇×Φ − Φ = (2) where k0 is the free-space wavenumber, and

Figure 1: Transverse cross-section of 3-D optical waveguide surrounded by PMLs.

[ ]

1 0 0 0 0 0 0 0 0 0 0 0 0 rx y z x x y ry z x y z rz x y z s s s p p p s s s p s s s µ µ µ −  ⋅        =_{ }= ⋅      _⋅     (3)

[ ]

0 0 0 0 0 0 0 0 0 0 0 0 rx y z x x y ry z x y z rz x y z s s s q q q s s s q s s s ε ε ε  ⋅        = = ⋅      _⋅     (4) for Φ =E and

[ ]

1 0 0 0 0 0 0 0 0 0 0 0 0 rx y z x x y ry z x y z rz x y z s s s p p p s s s p s s s ε ε ε −  ⋅        = = ⋅      _⋅     (5)

[ ]

0 0 0 0 0 0 0 0 0 0 0 0 rx y z x x y ry z x y z rz x y z s s s q q q s s s q s s s µ µ µ  ⋅        =_{ }= ⋅      _⋅   _ _ (6)

for Φ =H, where sx, sy, and sz are PML parameters as listed in Table 1 and the complex values of sj (j=1–4) are taken

as

1

j j

s = − jα . (7) The value of αj is used to control the attenuation of the electromagnetic field in PML regions and we choose it as

2 ,max j j j d ρ α =α _ _   (8)

where ρ is the distance from the beginning of the PML, dj (j=1–4) is the thickness of the PML as defined in Fig. 1, and

the subscript “max” refers to the maximum value.

To obtain the modal effective index, neff, and the mode field distribution, the field is split into a transverse and a

longitudinal component in the form

(

x y z, ,

)

φ_t

( )

x y, zφ_z

( )

x y, exp( j zβ )

Φ =_ +
_ − (9) where β = k0neff is the propagation constant. Dividing the waveguide cross-section into curvilinear hybrid edge/nodal

elements11 as shown in Fig. 2, the field components φx, φy, and φz in each element are expanded as

Table 1 PML parameters PML region PML parameter 1 2 3 4 5 6 7 8 sx s1 s2 1 1 s1 s2 s1 s2 sy 1 1 s3 s4 s3 s3 s4 s4 sz 1 1 1 1 1 1 1 1

Figure 2: Curvilinear hybrid edge/nodal element based on linear tangential and quadratic normal vector basis functions.

{ } { }

{ } { }

{ } { }

T t e x T y _{t e} T z z e U V j N φ φ φ φ φ φ _φ           =_{  }= _       _ _ (10)

where {φt}e and {φz}e are edge and nodal variables for each element, respectively, {U} and {V} are the shape function

vectors for edge elements, and {N} is the shape function vector for nodal elements.11

Using (9) and (10) in (2) and applying the finite element technique to the transverse x-y plane, we obtain the following matrix eigenvalue equation

[ ]

{ }

2 2

[ ]

{ } { }

0 eff 0 K φ −k n M φ = (11) with

{ } { }

_{{ }}

t z φ φ φ   =_{ }     (12)

[ ] [ ] [ ]

_{[ ] [ ]}

0 0 0 tt K K =_ _    ,

[ ] [ ] [ ]

[ ] [

]

tt tz zt zz M M M M M   =      (13)

[ ]

{ } { }

{ } { }

{ } { }

{ } { }

{ }{ }

{ }{ }

2 2 0 0 T T T T tt z z z z e T T x y U V V U U U V V K p p p p y x x y y y x x k q U U k q V V dxdy  _{∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂}  = + − −  ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂   + + _

∑∫∫

₍₁₄₎

[ ]

{ }{ }

T

{ }{ }

T tt y x e M = _p U U +p V V _dxdy  

∑∫∫

(15)

[ ]

tz y

{ } { }

T x

{ } { }

T e N N M p U p V dxdy x y  _{∂ ∂}    = +  ∂ ∂   

∑∫∫

(16)

[ ]

{ }{ }

T

{ }{ }

T zt y x e N N M p U p V dxdy x y  ∂ ∂  =  _∂ + _∂     

∑∫∫

(17)

[

]

{ } { }

{ } { }

2

{ }{ }

0 T T T zz y x z e N N N N M p p k q N N dxdy x x y y  _{∂ ∂ ∂ ∂}    = + −  ∂ ∂ ∂ ∂   

∑∫∫

(18)

where {0} is a null vector, [0] is a null matrix, and Σe extends over all different elements.

We solve the matrix eigenvalue equation (11) by the modified SIPM to obtain the eigenvalue, neff, and the

corresponding eigenvector, {φ}.

2.2. Finite element beam propagation method

Again referring to Fig. 1 and the related coordinate system, we now consider an anisotropic optical waveguide with the relative permittivity and the permeability tensors

[ ]

r xxyx xyyy xzyz zx zy zz ε ε ε ε ε ε ε ε ε ε     =        ,

[ ]

xx xy xz r yx yy yz zx zy zz µ µ µ µ µ µ µ µ µ µ     =        . (19)

Using PMLs for anisotropic media,10 the propagation of the field is described by the same vectorial wave equation (10) but with

[ ]

1 0 0 0 0 0 0 xx xy xz xx y z x yx yy yz yy z x y zx zy zz zz x y z p p p s s s p p p p s s s p p p s s s µ µ µ −    ⋅      =  = ⋅     _⋅          (20)

[ ]

xxyx xyyy xzyz xx z yxy z x yy z xyz x y x yzy xz zx zy zz y zx x zy zz x y z q q q s s s s s q q q q s s s s s q q q s s s s s ε ε ε ε ε ε ε ε ε    ⋅      =  = ⋅     _⋅          (21) for Φ =E and

[ ]

1 xx xy xz xx y z x z xy y xz yx yy yz z yx yy z x y x yz zx zy zz y zx x zy zz x y z p p p s s s s s p p p p s s s s s p p p s s s s s ε ε ε ε ε ε ε ε ε −    ⋅      =  = ⋅     _⋅          (22)

[ ]

0 0 0 0 0 0 xx xy xz xx y z x yx yy yz yy z x y zx zy zz zz x y z q q q s s s q q q q s s s q q q s s s µ µ µ    ⋅      =  = ⋅     _⋅          (23)

for Φ =H, where sx, sy, and sz are the same PML parameters as in Table 1.

Taking an appropriate reference refractive index n0 and assuming the slowly varying envelope approximation

(SVEA), we write

(

)

(

)

₍

₎

₍

₎

0 0 , , _x , , _y , , _z , ,
exp( ). x y z φ x y z x φ x y z y φ x y z z jk n z Φ =_ + + _ − (24)

Dividing the waveguide cross-section into curvilinear hybrid edge/nodal elements as shown in Fig. 2, the field components φx, φy, and φz in each element are again expanded as in (10).

Using (24) and (10) in (2), employing the longitudinal-field transformation12

(

_{, ,}

)

_{exp( 0 0 )}

{

₍

_{, ,}

₎

_{exp( 0 0 )}

}

z x y z jk n z j z x y z jk n z

z

φ − = ∂ φ −

∂ (25)

and applying the finite element technique to the transverse x-y plane, we obtain the following matrix equation:

[ ]

2

{ }

(

[ ]

[ ]

)

{ }

(

[ ]

[ ]

2 2

[ ]

)

{ } { }

0 0 0 0 0 0 2 2 0 d d M L jk n M K jk n L k n M dz dz φ φ φ + − + − − = (26) with

{ }

_{{ }}

{ }

t z φ φ φ     =     (27)

[ ] [ ] [ ]

_{[ ] [ ]}

0 0 0 tt K K =_ _    ,

[ ] [ ] [ ]

[ ] [ ]

0 tt tz zt L L L L   =_{ }    ,

[ ] [ ] [ ]

[ ] [

]

tt tz zt zz M M M M M   =_{ }     (28)

[ ]

{ } { }

{ } { }

{ }{ }

{ }{ }

{ } { }

{ } { }

_{{ }{ }}

_{{ }{ }}

2 2 0 0 2 2 0 0 T T T T tt zz zz xx xy e T T T T zz zz yx yy U V U U K p p k q U U k q U V y x y y V V V U p p k q V U k q V V dxdy x x x y  _{∂ ∂ ∂ ∂}  = − + +  ∂ ∂ ∂ ∂   ∂ ∂ ∂ ∂ _ − + + +  ∂ ∂ ∂ ∂ _

∑∫∫

(29)

[ ]

{ } { }

{ } { }

{ }{ }

{ }{ }

{ } { }

{ } { }

{ }{ }

{ }{ }

T T T T tt yz yz zx zy e T T T T xz xz zx zy V U U U L p U p U p V p U x y y y V U V V p V p V p V p U dxdy x y x x  _{∂ ∂ ∂ ∂}  = − − +  ∂ ∂ ∂ ∂   ∂ ∂ ∂ ∂ _ − + + −  ∂ ∂ ∂ ∂ _

∑∫∫

(30)

[ ]

{ } { }

{ } { }

{ } { }

{ } { }

{ }{ }

{ }{ }

2 2 0 0 T T T T tz zx zy zx zy e T T xz yz U N U N V N V N L p p p p y y y x x y x x k q U N k q V N dxdy  _{∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂}  = − + + −  ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂   − − _

∑∫∫

₍₃₁₎

[ ]

{ } { }

{ } { }

{ } { }

{ } { }

{ }{ }

{ }{ }

2 2 0 0 T T T T zt xz yz xz yz e T T zx zy N U N U N V N V L p p p p y y x y y x x x k q N U k q N V dxdy  _{∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂}  = − − +  ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂   + + _

∑∫∫

₍₃₂₎

[ ]

{ }{ }

T

{ }{ }

T

{ }{ }

T

{ }{ }

T tt yx yy xx xy e M = _−p U V +p U U +p V V −p V U _dxdy  

∑∫∫

(33)

[ ]

tz yx

{ } { }

T yy

{ } { }

T xx

{ } { }

T xy

{ } { }

T e N N N N M p U p U p V p V dxdy y x y x  _{∂ ∂ ∂ ∂}    = − + + −  ∂ ∂ ∂ ∂   

∑∫∫

(34)

[ ]

{ }{ }

T

{ }{ }

T

{ }{ }

T

{ }{ }

T zt xy yy xx yx e N N N N M p U p U p V p V dxdy y x y x  ∂ ∂ ∂ ∂  = − _∂ + _∂ + _∂ − _∂     

∑∫∫

(35)

[

]

{ } { }

{ } { }

{ } { }

{ } { }

{ }{ }

2 0 . T T T T zz yx yy xx xy e T zz N N N N N N N N M p p p p x y x x y y y x k q N N dxdy  _{∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂}  = − + + −  ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂   − _

∑∫∫

₍₃₆₎

Equation (26) is the basic equation for the beam propagation analysis. Under the Fresnel approximation, (26) becomes

[ ]

[ ]

(

)

{ }

(

[ ]

[ ]

2 2

[ ]

)

{ } { }

0 0 0 0 0 0 2 d 0 L jk n M K jk n L k n M dz φ φ − + − − = . (37)

For waveguides with diagonal permittivity and permeability tensors, εkl=0 and µkl=0 (k≠l), and (37) reduces to

[ ]

{ }

(

[ ]

2 2

[ ]

)

{ } { }

0 0 0 0 2jk n M d K k n M 0 dz φ _φ − + − = . (38) Applying the Crank-Nicholson algorithm, we obtain

[ ]

Ai

{ }

φ i+1=

[ ]

B i

{ }

φ i (39)

with

[ ]

A_i = −2jk n_{0 0}

[ ]

M +0.5∆z

(

[ ]

K −k n_{0 0}2 2

[ ]

M

)

(40) and

[ ]

B _i = −2jk n_{0 0}

[ ]

M −0.5∆z

(

[ ]

K −k n_{0 0}2 2

[ ]

M

)

. (41) 2.3. Imaginary-distance propagation method

Our purpose is to seek the eigenmodes of the optical waveguides from the propagation scheme shown above and an initial field is needed to start the process. Assuming m eigenmodes including radiation modes propagate in the waveguides, at the ith propagation step the field {φ}i can be expressed as a superposition of the eigenmodes

{ }

,

{ }

1 m r i r i r A f φ = =

∑

(42) where {fr} is the field distribution of the rth eigenmode and Ar,i indicates the complex amplitude of the rth eigenmode.

On the other hand, from the earlier discussion of the basics of the eigenmode solver, we know that the effective refractive index, neff,r, and the corresponding field, {fr}, of the rth eigenmode must satisfy the following matrix

eigenvalue equation

[ ]

{ }

2 2

[ ]

{ } { }

0 , 0

r eff r r

K f −k n M f = . (43) Substituting (43) in (39)–(41), the field distribution of the rth eigenmode after a propagation distance of ∆z

becomes

{ }

(

₍

)

₎

{ }

2 2 2 0 0 0 , 0 1 2 2 2 0 0 0 , 0 2 0.5 2 0.5 eff r r _i r _i eff r jk n zk n n f f jk n zk n n + − − ∆ − = − + ∆ − . (44)

Observing the equation above, the propagation step size along the z-direction, ∆z, can be suitably chosen as

(

2 0 2

)

, 0 0 4
eff r n z j n n k ∆ − (45) so that the rth eigenmode can see a very large gain and after a sufficiently large number of i steps, {φ} will converge to the rth eigenvector {fr},13 The effective index of this eigenmode, neff,r, is obtained by

{ }

[ ]

{ }

{ }

[ ]

{ }

† 2 , _{2 †} 0 i i i eff ri i i i K n k M φ φ φ φ = (46) at the ith propagation step.

3. NUMERICAL RESULTS

We first examine if the ID-BPM and the eigenmode solver give the same results. Fig. 3(a) shows the cross-section of a microstructured fiber, or a holey fiber, with three rings of 36 air holes in a hexagonal lattice, as studied in Ref. 5, where d is the air hole diameter and the lattice pitch Λ is taken to be 2.3 µm. Fig. 3(b) shows the confinement loss in dB/m, which is obtained from the imaginary part of the modal propagation constant, of the x-polarized fundamental mode as a function of the wavelength for three different d/Λ values. It is seen that very good agreement exists among our FE-ID-BPM and modified SIPM (MSIPM) (eigenmode solver) results and the FE-ID-FE-ID-BPM results of Ref. 5. Figure 4(a) and (b) shows the corresponding results for the cases with two rings of 18 air holes and one ring of 6 air holes, respectively. Again very good agreement exists among these methods. We then consider the PBG fiber with honeycomb shape analyzed in Ref. 6, as shown in Fig. 5(a), where we depict one quarter of a fiber of nine hole rings, whereas the fiber of two hole rings is completely shown. There is an extra hole at the center of the fiber, representing a defect. All the holes have the same diameter d. We take Λ = 1.62 µm and the wavelength λ = 1.55 µm as in Ref. 6. The element division profile of the nine-ring case is shown in Fig. 5(b). For a profile with so many elements used, the MSIPM (eigenmode solver) needs much time to obtain the results, and we employ our FE-ID-BPM to analyze this structure. Figure 5(c) shows the comparison of our FE-ID-BPM calculation of the confinement loss versus the wavelength with that of Ref. 6 for d/Λ = 0.41 and 0.55 and the number of hole rings going from 1 to 9. Good agreement is observed for the structures with the number of hole rings less than 4. However, our calculation becomes obviously different from that of Ref. 6 for the number of hole rings larger than 6. It is seen that our results vary more smoothly with the number of rings, while those of Ref. 6 show a sudden change in the decreasing rate near the location where the number of rings is six. We thus believe our results are more reasonable.

(a) (b) Figure 3: (a) Cross-section of a microstructured fiber with three rings of 36 air holes. (b) Wavelength

dependence of the confinement loss of the x-polarized fundamental mode in the fiber of (a).

(a) (b) Figure 4: (a) Same as Fig. 3(b) but with two rings of 18 air holes. (b) Same as Fig. 3(b) but with one

ring of 6 air holes.

(a) (b) (c)

Figure 5: (a) Cross-section of honeycomb PBG fibers. (b) Element division profile of a quarter of a honeycomb fiber with nine rings. (c) Confinement loss versus the number of hole rings for the fiber of (a) at λ = 1.55 µm for d/Λ = 0.41 and 0.55.

4. CONCLUSION

We have developed a full-vectorial imaginary-distance BPM based on curvilinear hybrid edge/nodal elements with triangular shape as in Ref. 5 and a related eigenmode solver for the analysis of microstructured and PBG fibers, including the determination of the modal confinement losses. A modified shift inverse power method is used for solving the eigenvalue problem. The two solution schemes are found to give consistent results. Our calculation shows some

difference in the confinement loss values from those reported in Ref. 6 for the honeycomb PBG fibers. From the behavior of the loss curve versus the number of hole rings, we believe our results are more reasonable.

Acknowledgements: This work was supported in part by the National Science Council of the Republic of China under grants NSC92-2215-E-002-008 and NSC93-2215-E-002-042.

5. REFERENCES

1. J. C. Knight, T. A. Birks, P. St. J. Russell, and D. M. Atkin, “All-silica single-mode optical fiber with photonic crystal cladding,” Opt. Lett., vol. 21, pp. 1547-1549, 1996.

2. P. St. J. Russell, “Photonic crystal fibers,” Science, vol. 299, pp. 358-362, 2003.

3. J. C. Knight, J. Broeng, T. A. Birks, and P. St. J. Russell, “Photonic band gap guidance in optical fibers,” Science, vol. 282, pp. 1476-1478, 1998.

4. T. P. White, R. C. McPhedran, C. M. de Sterke, L. C. Botton, M. J. Steel, “Confinement losses in microstructured optical fibers,” Opt. Lett., vol. 26, pp. 1660-1662, 2001.

5. K. Saitoh, and M. Koshiba, "Full-vectorial imaginary-distance beam propagation method based on a finite element scheme: application to photonic crystal fibers," J. Quantum Electronics, vol. 38, pp. 927-933, 2002.

6. D. Ferrarini, L. Vincetti, and M. Zobobi, "Leakage properties of photonic crystal fibers," Opt. Express, vol. 10, pp. 1314-1319, 2002

7. D. Yevick and B. Hermansson, “New approach to lossy optical waveguide,” Electron. Lett., vol. 21, pp. 1029-1030, 1985.

8. S. Selleri, L. Vincetti, A. Cucinotta, M. Zoboli, "Complex FEM modal solver of optical waveguides with PML boundary conditions," Opt. Quantum Electron., vol. 33, pp. 359-371, 2001.

9. S. Rebay, “Efficient unstructured mesh generation by means of Delaunay triangulation and Bowyer-Watson algorithm,” J. Comput. Phys., vol. 106, pp. 125-138, 1993.

10. F. L. Teixeira, and W. C. Chew, "General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media," IEEE Microwave Guided Wave Lett., vol. 8, pp. 223-225, 1998.

11. M. Koshiba, and Y. Tsuli, "Curvilinear hybrid edge/nodal elements with triangular shape for guided-wave problems," J. Lightwave Technol., vol. 18, pp. 737-743, 2000.

12. D. Schulz, C. Gingener, M. Bludsuweit, and E. Voges, "Mixed finite element beam propagation method," J.

Lightwave Technol., vol. 16, pp. 1336-1341, 1998.

13. C. L. Xu, W. P. Huang, and S. K. Chaudhuri, "Efficient and accurate vector mode calculations by beam propagation method," J. Lightwave Technol., vol. 11, pp. 1209-1215, 1993.