適於光電研究的先進有限差分相關數值方法與模型之發展(2/3)

行政院國家科學委員會專題研究計畫 期中進度報告

適於光電研究的先進有限差分相關數值方法與模型之發展

(2/3)

期中進度報告(精簡版)

計 畫 類 別 ： 個別型

計 畫 編 號 ： NSC 95-2221-E-002-380-

執 行 期 間 ： 95 年 08 月 01 日至 96 年 07 月 31 日

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

計 畫 主 持 人 ： 張宏鈞

報 告 附 件 ： 出席國際會議研究心得報告及發表論文

處 理 方 式 ： 期中報告不提供公開查詢

中 華 民 國 96 年 05 月 31 日

行政院國家科學委員會專題研究計畫期中進度報告

適於光電研究的先進有限差分相關數值方法與模型之發展(2/3

)

Development of Advanced Finite Difference Based Numerical

Methods and Models for Photonics Research (2/3)

計畫編號：NSC 95-2221-E-002-380

執行期限：95 年 8 月 1 日至 96 年 7 月 31 日

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

計畫參與人員：江柏叡(台大光電所) 陳明昀(台大光電所)

林彥宏(台大光電所) 張峰偉(台大電信所)

江舜凡(台大光電所)

摘要

本計畫發展數種有限差分相關的頻域與時域數值模型，以做為研究各種光導波元件

的工具。

第二年至目前的執行情形，依

有限差分時域法模型之發展、以有限差分頻域法計

算三維光子晶體與各向異性光子晶體能隙結構、譜方法模型之發展與應用、包含光子晶

體光纖之特殊光波導之分析與研究等四項研究課題，分述於本報告。

關鍵詞：光波導、光子晶體、光子晶體波導、波導元件、有限差分頻域法、有限差分時

域法、時域譜方法、頻域譜方法

Abstract

This research concerns development of several finite difference (FD) based

frequency-domain and time-domain numerical models as tools for studying various photonic

guided-wave devices. In this report we describe the research results up to date in the second

year of this three-year project, including the development of the finite-difference time-domain

(FDTD) method, the analysis of 3-D photonic crystal and anisotropic photonic crystal band

structures using the finite-difference frequency-domain (FDFD) method, the development and

applications of pseudo-spectral method models, and the analysis and investigation of photonic

crystal fibers and special-type optical waveguides.

Keywords: Optical waveguides, photonic crystals, photonic crystal waveguides, optical

waveguide devices, finite-difference frequency-domain method, finite-difference

time-domain method, pseudo-spectral time-domain method, pseudo-spectral

frequency method

I. Development of the FDTD Method

We have established both two-dimensional (2-D) and 3-D FDTD numerical models for

studying various photonic structures. For 3-D FDTD modeling, the simulation is quite CPU

time consuming and memory demanding and going for parallel computing is a natural way.

In the first year, we had set up a personal computer cluster in our laboratory for such

simulation and had made our 3-D FDTD code a parallel computing model [1]. We have

continuously extended the capacity of our models to handle more general and more

complicated structures, for example, the inclusion of material dispersion models and the

treatment of nonlinear optics problems. We briefly show two numerical examples in the

following.

The first example relates to the subwavelength imaging using an array of silver nanorods,

as reported in [2]. The left panel of Fig. 1 is Fig. 1 of [2], which shows the nanorod-array

structure and the FDTD simulation results of the image profiles at different places due to

point sources shaped as the letter “λ” as an object. The rod height is 50 nm, the rod diameter

is 20 nm, and the side-to-side distance between the rods is 40 nm. The right panel of Fig. 1

shows our 3-D FDTD simulation results corresponding to (f) and (g) places in the left panel.

Our FDTD is with the Drude dispersion model and good agreement with {2} is seen. Fig. 2(a)

is from Fig. 2 of [2], showing the light intensity distribution along one rod, and the blue line

in Fig. 2(b) is our calculation. Again, excellent agreement with [2] is observed.

(Left panel) (Right panel)

Fig. 1 Left panel: Fig. 1 of [2]. Right panel: simulation results using our 3-D FDTD model.

(a) (b)

Fig. 2 (a) From Fig. 2 of [2]: light intensity distribution along one rod. (b) The

blue line is from our

FDTD calculation.

The second example relates to the FDTD simulation of nonlinear medium phenomena.

We examine our established model using a four-wave mixing problem in the same way as

was recently done in [3]. The simulation arrangement is as shown in Fig. 3(a), where uniform

plane pump and signal waves are incident from a linear medium upon a nonlinear medium

with 3

rd

-order instantaneous Kerr nonlinearity represented by the nonlinear susceptibility χ

0(3)

= 1.0 x 10

-18

m

2

/V

2

, The linear index of both media is 1.2. The pump wave is at 192 THz and

the signal wave is at 195 THZ, with their electric field amplitudes being 5.1 x 10

7

V/m and

5.1 x 10

6

V/m, respectively. The output spectrum of the electric field obtained from our

FDTD simulation is given in Fig. 3(b), showing converted waves at 189 and 198 THz, which

is consistent with the results of [3] and the analytical four-wave mixing theory.

(a) (b)

Fig. 3 (a) FDTD simulation arrangement for a four-wave mixing problem. (b)

The output spectrum

of the electric field obtained from our FDTD simulation.

II. Analysis of 3-D Photonic Crystal and Anisotropic Photonic Crystal Band Structures

Using the FDFD Method

We have extended our FDFD photonic crystal solver based on the Yee mesh [4] to the

analysis of 3-D photonic crystals (PCs) as well as PCs involving anisotropic materials. We

have demonstrated that the FDFD method can be efficiently utilized to calculate the band

diagrams of 3-D PC structures [5], [6]. We have considered two structures: the dielectric

spheres in air and the scaffold structure [5]. When the index averaging scheme is used, the

obtained band diagrams agree well with those calculated by the plane-wave-based transfer

matrix method and the plane-wave expansion method. The established FDFD method can

treat PCs with diagonal permittivity tensor. Please refer to Appendix I [5] for the detail. We

have also developed an Yee-mesh-based FDFD method for calculating the band structures of

2-D PCs containing anisotropic materials [7]. As a numerical example, we have analyzed a

square-lattice PC composed of circular cylinders of nematic liquid-crystal (LC) material with

silicon as background. The LC director is assumed to be perpendicular to the PC extension

direction so that the E and H waves are decoupled. The tenability of the band-gap size has

been studied and symmetry properties of the gap size versus the orientation of the LC director

observed [7]. Please refer to Appendix II [7].

III. Development and applications of pseudo-spectral method models

We have been working on the development of electromagnetic numerical models based

on pseudo-spectral methods both in frequency domain and in time domain. A novel analysis

method based on a multidomain pseudospectral method has been proposed for calculating the

band diagrams of 2-D PCs and is shown to possess excellent numerical convergence behavior

and accuracy [8]. The proposed method shows uniformly excellent convergence

characteristics for both the transverse-electric and transverse-magnetic waves in the analysis

of different structures. The analysis of a mini band gap is also shown to demonstrate the

extremely high accuracy of the proposed method. The first page of the paper published in

Physical Review E [8] is attached as Appendix III.

In the development of the pseudospectral time-domain (PSTD) simulation method, the

Drude dispersion model has been incorporated into the code for treating problems involving

metallic materials in the optical frequency domain. In related work, optical properties of a

metallic sphere based on the classical Drude model have been studied theoretically. The

comparison between the metallic sphere and the perfectly conducting one shows conspicuous

distinction in the electric field amplitude on the surface of the sphere and in the radar cross

section behavior [9]. Please refer to Appendix IV [9].

IV. Analysis and Investigation of Photonic Crystal Fibers and Special-type Optical

Waveguides

The FDFD waveguide eigenmode solver based on the Yee mesh [4] is a versatile

method for the analysis of different optical waveguides. We have used it to study the

propagation characteristics and chromatic dispersion coefficients of photonic crystal fibers

[10]–[12]. We have recently applied the FDFD solver to investigate modal properties of

guided modes on air-core terahertz waveguides with the air core surrounded by Teflon rods.

The antiresonant reflection optical waveguide (ARROW) model was successfully adopted to

explain the guiding behavior [13]. Please refer to Appendix V [13]. Based on such FDFD

analysis, we have proposed a very simple procedure for numerical calculation of chromatic

dispersion coefficients of photonic crystal fibers, which involves Chebyshev-Lagrange

interpolation polynomials [14]. No direct numerical differentiation of the effective refractive

index is needed. Please refer to Appendix VI [14].

References

[1] M. F. Chen and H. C. Chang, “Finite-Difference Time-Domain Analysis of Three-Dimensional Photonic-Wire Bends via Parallel Computing,” in Proceedings of Optics and Photonics Taiwan ’05 (OPT ’05) (CD-ROM), paper B-SA-IV10-5, Tainan, Taiwan, R.O.C., December 9–10, 2005.

[2] [3]

[4]

A. Ono, J. Kato, and S. Kawata, “Subwavelength optical imaging through a metallic nanorod array,” Phys. Rev. Lett., vol. 95, 267407, 2005.

M. Fujii, C. Koos, C. Poulton, I. Sakagami, J. Leuthold, and W. Freude, “A simple and rigorous verification technique for nonlinear FDTD algorithms by optical parametric four-wave mixing,” Microwave and Optical Technology Letters, vol. 48, pp. 88–91, 2006.

C. P. Yu and H. C. Chang, “Chapter 7: Applications of the Finite Difference Frequency Domain Mode Solution Method to Photonic Crystal Structures (50 pages),” in Electromagnetic Theory and Applications for Photonic Crystals (448 pages), edited by Kiyotoshi Yasumoto, pp. 351–400. Boca Raton, Florida: Marcel Dekker/CRC Press, Inc., 2006. (ISBN: 0849336775)

[5] Y. H. Lin, C. P. Yu, and H. C. Chang, “Calculation of Three-Dimensional Photonic-Crystal Band Diagrams Using the Finite-Difference Frequency-Domain Method,” in Optics and Photonics Taiwan ’06 (OPT ’06) (CD-ROM), paper AO-61, Hsinchu, Taiwan, R.O.C., December 15–16, 2006.

H. C. Chang, Y. H. Lin, C. P. Yu, and M. M. Chen, “Finite Difference Analysis of 3D Anisotropic Photonic Crystals, in Technical Digest, PECS-VII: International Symposium on Photonic and Electromagnetic Crystal Structures VII (CD ROM), paper A21, Monterey, California, April 8–11, 2007. [6]

M. M. Chen, S. M. Hsu, and H. C. Chang, “Analysis of 2-D Photonic Crystals Involving Liquid Crystals Using the Finite-Difference Frequency-Domain Method,” in OSA 2006 Integrated Photonics and [7]

Nanophotonics Research and Applications (IPNRA ’07) Technical Digest (CD ROM), paper IMB3 (3 pages), Salt Lake City, Utah, July 8–11, 2007.

P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of Two-Dimensional Photonic Crystals Using a Multidomain Pseudospectral Method,” Physical Review E, Vol. 75, pp. 026703-1–026703-14, 20 February 2007.

[8]

B. Y. Lin, C. H. Teng, and H. C. Chang, “Near-Field and Far-Field Behavior of a Metallic Nanoscale Sphere at Optical Frequencies Based on the Classical Drude Model,” accepted for presentation at URSI Commission B – International Symposium on Electromagnetic Theory (EMTS 2007), Ottawa, ON, Canada, August 26–28, 2007.

[9]

[10] H. C. Chang, C. P. Yu, P. J. Chiang, and S. M. Hsu, “Development and Applications of High-Accuracy Optical Waveguide Eigenmode Solvers,” presented at International Workshop on Photonics and Display Technologies, National Taiwan University of Science and Technology, Taipei, Taiwan, R.O.C., September 22–23, 2006. (Invited paper)

[11] H. C. Chang, “High-Accuracy Analysis of New-Type Optical Waveguide Structures,” Optics and Photonics Taiwan ’06, Hsinchu, Taiwan, R.O.C., December 15–16, 2006. (Invited paper)

[12] H. C. Chang and C. P. Yu, “Finite-Difference Numerical Analysis of Photonic Crystals and Photonic Crystal Fibers,” presented at Physical Society R.O.C. 2007 Annual Meeting, ChungLi, Taiwan, R.O.C., January 23–25, 2007. (Invited paper)

[13] C. P. Yu and H. C. Chang, “Air-Core Waveguides for Terahertz Signals,” in OSA 2006 Integrated Photonics and Nanophotonics Research and Applications (IPNRA ’07) Technical Digest (CD ROM), paper ITuH4 (3 pages), Salt Lake City, Utah, July 8–11, 2007.

[14] P. J. Chiang, C. P. Yu, and H. C. Chang, “Robust Calculation of Chromatic dispersion Coefficients of Optical Fibers from Numerically Determined Effective Indices Using Chebyshev-Lagrange Interpolation Polynomials,” IEEE/OSA Journal of Lightwave Technology, Vol. 24, No. 11, pp. 4411–4416, November 2006.

CALCULATION OF THREE-DIMENSIONAL

PHOTONIC-CRYSTAL BAND DIAGRAMS

USING THE FINITE-DIFFERENCE FREQUENCY-DOMAIN METHOD

Yen-Hung Lin (林彥宏), Chin-Ping Yu (于欽平),and Hung-Chun Chang* (張宏鈞)

Graduate Institute of Electro-Optical Engineering, National Taiwan University,

Taipei, Taiwan 106-17, R.O.C.

*also with Department of Electrical Engineering and

Graduate Institute of Communication Engineering, National Taiwan University

Phone: +886-2-23635251-513, Fax: +886-2-23638247, E-mail: [email protected]

(NSC94-2215-E-002-022 and NSC95-2221-E-002-380)

Abstract ---

We adopt the finite-difference frequency-domain method to determine the 3D photonic crystal band structures. The band diagrams of two different kinds of structures are efficiently obtained with good agreement with results from other methods.

Keywords:

Finite-difference frequency-domain method, photonic band gap, photonic crystals.

INTRODUCTION

Photonic crystals (PCs) have many applications, including the routing and guiding of light in all three dimensions, due to the photonic band gaps (PBGs) resulting from periodic structures [1]. A simple cubic (sc) lattice structure is symmetrical with respect to x, y, and z directions and therefore is the most natural platform for constructing photonic crystal devices [2]. Besides, theoretical investigation suggests that it is possible to create such an sc structure with a complete PBG [3]. In this paper, we develop a finite-difference frequency-domain (FDFD) method based on the Yee mesh to calculate the band diagrams of three-dimensional (3D) PCs. Two-dimensional version of the FDFD method has been used to calculate the band structures of 2D PCs [4]. We will consider two kinds of 3D structures. The first one is dielectric spheres in air, and the second is scaffold structures, as shown in Fig. 1(a) and (b), respectively. The band diagrams of these two structures are obtained and compared with other numerical methods and experimental results.

Hx( i +0.5, j,k) Hz( i ,j,k+0.5) Ey( i +0.5,j,k+0.5) ( i ,j,k) Hz ( i +1,j ,k+0.5) Hx(m+0.5,n,k+1) Ez( i +0.5,j +0.5,k) Ez( i +0.5,j+0.5,k+1) Hy( i +1,j +0.5,k) Hy( i, j+0.5,k) Hy( i ,j+0.5,k+1) Hy( i +1,j+0.5,k+1) Ex( i, j+0.5,k+0.5) Ex( i +1,j +0.5,k+0.5) Ey( i +0.5,j +1,k+0.5) Hz( i +1,j+1,k+0.5) Hz( i ,j +1,k+0.5) Hx( i +0.5,j +1,k) Hx( i +0.5,j +1,k+1) (a) (b) (c)

Figure 1. (a) The 3D sc photonic crystal structure composed of dielectric spheres in air. (b) The 3D sc photonic crystal structure composed of scaffold structures. (c) 3D Yee’s mesh for the FDFD method.

NUMERICAL METHOD

Since 3D PCs are periodic along three orthogonal directions, the propagation in all directions must be taken into account, i.e., the three phase constants kx, ky, and kz may all have non-zero values. Starting from

Maxwell’s curl equations and using the central difference scheme with 3D Yee’s mesh shown in Fig. 1(c), we can obtain the following six equations with six electromagnetic components as

(3) (2) (1) ) , , ( , ) , , ( , ) , 1 , ( , ) , 1 , ( , ) , 1 , ( , 0 ) 1 , , ( , ) 1 , , ( , ) , , ( , ) , , ( , ) 1 , , ( , 0 ) , 1 , ( , ) , 1 , ( , ) 1 , , ( , ) 1 , , ( , ) 1 , 1 , ( , 0 2 1 2 1 2 1 2 3 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 3 2 1 2 1 2 1 2 1 2 3 2 1 2 1 2 1 2 3 2 1 2 1 y E E x E E H j x E E z E E H j z E E y E E H j k j i x k j i x k j i y k j i y k j i z k j i z k j i z k j i x k j i x k j i y k j i y k j i y k j i z k j i z k j i x ∆ − − ∆ − = − ∆ − − ∆ − = − ∆ − − ∆ − = − + + + + + + − + + + + + + + − + + + + + + + + + + + − + + − + + − + + − + + − ωµ ωµ ωµ z y x

Appendix I

) 6 ( ) 5 ( ) 4 ( ) 1 , , ( , ) 1 , 1 , ( , ) 1 , , 1 ( , ) 1 , , ( , ) 1 , , ( , 0 ) , 1 , 1 ( , ) , 1 , ( , ) , 1 , ( , ) 1 , 1 , ( , ) , 1 , ( , 0 ) , , ( , ) 1 , , ( , ) , , ( , ) , 1 , ( , ) , , ( , 0 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 y H H x H H E j x H H z H H E j z H H y H H E j k j i x k j i x k j i y k j i y k j i z z k j i z k j i z k j i x k j i x k j i y y k j i y k j i y k j i z k j i z k j i x x ∆ − − ∆ − = ∆ − − ∆ − = ∆ − − ∆ − = + − + + − + + − + + + + − + + − + + + − + + − + + − + + + + + + + + ε ωε ε ωε ε ωε

where ω is the angular frequency, ε0 and µ0 are the permittivity and permeability of free space, respectively, and

εx, εy, and εz are the relative permittivities at the corresponding grid points. Here, we have assumed the material

of the PC has permeability µ0, but can have a diagonal permittivity tensor. Therefore, more general anisotropic

materials can be considered by equations (1)−(6). Those equations can then be written into the matrix form. After some mathematical work, we can obtain an eigenvalue equation as

) 7 ( 1 1 1 2 0               =                             − − − − − −               z y x z y x x x y y z y z x y z x x z z y x x z x y y y z z zz yy xx E E E E E E ) U V U V ( ) U (V ) U (V ) U (V ) U V U V ( ) U (V ) U (V ) U (V ) U V U V ( ε 0 0 0 ε 0 0 0 ε k

where k02 = ω2µ0ε0, 0 is the zero matrix, εxx, εyy, and εzz are diagonal matrices representing εx, εy, and εz at the

corresponding grid points, and Ux, Uy, Uz, Vx, Vy, and Vz are square matrices determined by (1)−(6) and the

periodic boundary conditions at the edge of the computing window. By finding out the eigen frequencies of (7), we can construct the band diagrams of the PC structures.

NUMERICAL PROCESS AND RESULTS

The first 3D PC under study is the sc structure containing dielectric spheres in air. The sphere has a refractive index of n = 3.4 and a radius of r = 0.2a with a being the lattice constant of the PC. The unit cell of the structure is shown in Fig. 2(a). The reciprocal lattice of this PC and some high-symmetry points are schematically displayed as the inset of Fig. 2(b). The band diagram is plotted along the high-symmetry line Γ-X-M-R. For this simple cubic structure, we have carried out the FDFD calculation to get the eigen frequencies of (7). Numerical computation of the photonic band structure uses 16×16×16 grid points. In order to solve the discontinuity at dielectric interfaces, the relative permittivity ε(i, j, k) at each mesh point is determined by the index averaging scheme through ε(i, j, k) = fa εa + (1－fa) εd where fa is the filling fraction of air in the mesh, and

εa and εd are the relative permittivities of air and the dielectric medium, respectively. The band diagram of this

structure is shown in Fig. 2(b) with the dashed lines being our results from the FDFD method and the dots being from the plane-wave-based transfer-matrix method (TMM) [5]. It can be seen that our results agree quite well with those of [5]. However, no complete band gap can be found in this structure.

a r r=0.3a n=3.4 0 0.1 0.2 0.3 0.4 0.5 0.6 Γ X M R N orm al iz ed fre que nc y ( ω a/2 πc) (a) (b)

Figure 2. (a) The unit cell of the 3D sc PC containing dielectric spheres in the air. (b) The band diagram of the PC of (a) with the dashed lines obtained by the FDFD method and the dots representing the results from the plane-wave-based TMM [5]. Plane-wave-based TMM [5] FDFD method kx ky kz kx ky kz X R M Γ

We proceed further to consider a more complex 3D PC structure: the scaffold PC which consists of layers of interconnected square rods joined by vertical posts to form a sc lattice. The reciprocal lattice of the scaffold structure is the same as the inset of Fig 2(b). The unit cell of the scaffold structure is shown in Fig. 3(a), and the width of the square rod is w = 0.13950a. The dielectric filling fraction is 19% and the refractive index of dielectric materials is n = 3.6. By applying the FDFD method with 16×16×16 grid points in the unit cell and the index averaging scheme for the permittivity at each grid point, the band diagram of this structure can then be obtained as the dashed lines as shown in Fig. 3(b). The triangular dots in Fig. 3(b) are from the plane-wave-expansion (PWE) method [6], and the dots are from experimental results [2].

a w=0.13950a n=3.6 w 0 0.1 0.2 0.3 0.4 0.5 Γ X M R Experimental results [2] FDFD method (this work) PWE method [6] N orm al iz ed fr eq ue nc y ( ω a/2 πc) (a) (b)

Figure 3. (a) The unit cell of the 3D sc PC containing the scaffold structure. (b) The band diagram of the PC of (a) with dashed lines from the FDFD method, triangular dots from the PWE method [6], and the dots from experimental results [2].

The lowest six bands are presented in Fig. 3(b) with bands 1 and 2 being the photonic valence bands (VBs), and bands 3 to 6 being the photonic conduction bands (CBs). One can see a complete band gap in Fig. 3(b), which is defined by the VB maximum at the R point and the CB minimum at the X point. The lower bound of the band gap is ω1 = 0.367, and the upper bound is ω2 = 0.390. In our numerical calculation, there is a visible

difference between our results and those of [6] if the discontinuity of dielectric interfaces is not taken into account. By adopting the index averaging scheme, our results match well with those obtained by the PWE method [6], as shown in Fig. 3(b). Besides, from Fig. 3(b), it is seen that both our FDFD results and those of the PWE method are close to the trend of the dots from experimental results [2] at the X-M region in the VB region.

CONCLUSION

We have demonstrated that the FDFD method can be efficiently utilized to calculate the band diagrams of 3D PC structures. We have considered two structures: the dielectric spheres in air and the scaffold structure. When the index averaging scheme is used, the obtained band diagrams agree well with those calculated by the plane-wave-based TMM and the PWE method. The established FDFD method can treat PCs with diagonal permittivity tensor.

REFERENCES

[1] E. Yablonovitch, "Inhibited spontaneous emission in solid-state physics and electronics," Phys. Rev. Lett.

58, pp. 2059-2062 (1987).

[2] S. Y. Lin, J. G. Fleming, and R. Lin, "Complete three-dimensional photonic bandgap in a simple cubic structure," J. Opt. Soc. Am. B 18, pp. 32-35 (2001).

[3] H. S. Sozuer and J. W. Haus, "Photonic bands: simple cubic lattice," J. Opt. Soc. Am. B 10, pp. 296-302 (1993).

[4] C. P. Yu and H. C. Chang, "Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals," Optics Express 12, pp. 1397-1408 (2004)

[5] Z. Y. Li and L. L. Lin, "Photonic band structures solved by a plane-wave-based transfer-matrix method," Phys. Rev. E 67, 046607 (2003)

[6] J. Mizuguchi, Y. Tanaka, S. Tamura, and M. Notomi, "Focusing of light in a three-dimensional cubic photonic crystal," Phys. Rev. B 67, 075109 (2003)

ω2

Analysis of 2-D Photonic Crystals Involving Liquid

Crystals Using the Finite-Difference

Frequency-Domain Method

Ming-mung Chen, Sen-ming Hsu, and Hung-chun Chang∗

Graduate Institute of Electro-Optical Engineering National Taiwan University, Taipei, Taiwan 106-17, R.O.C.

∗_{also with Department of Electrical Engineering and Graduate Institute of}

Communication Engineering, National Taiwan University

Phone: +886-2-23635251-513 Fax: +886-2-23683824 E-mail: [email protected]

Abstract: The finite-difference frequency-domain method is formulated for calculating band structures of 2-D photonic crystals involving anisotropic materials such as liquid crys-tals. The director of the liquid crystal is allowed to rotate in the plane of the unit cell.

c

°2007 Optical Society of America

OCIS codes: (160.3710) Liquid crystals; (260.2110) Electromagnetic theory; (999.9999) Pho-tonic crystals.

1. Introduction

Research on photonic crystals (PCs) has been actively conducted for two decades since the pioneering works by Yablonovitch [1] and John [2] regarding that a periodic dielectric structure exhibits a forbidden band for optical energy. Many applications of PCs have been demonstrated, including novel waveguiding structures, waveguide devices, cavities, filters, multiplexers, etc. Several studies have been performed on PCs made of anisotropic materials. Zabel and Stroud [3] first analyzed three-dimensional (3-D) anisotropic PCs and showed the narrowing of the band gap. Li et al. [4] demonstrated that for 2-D PCs, the complete band gap can be achieved by including uniaxial materials with the optical axis along the extension direction of the PC. In this paper, we particularly consider the more involved structure with the optical axis perpendicular to the PC extension direction.

The analysis of PCs and their photonic band gaps (PBGs) have been most often based on the plane-wave expansion (PWE) method, including the study of anisotropic PCs. In this research, we develop a finite-difference frequency-domain (FDFD) method based on the Yee mesh for calculation of band struc-tures of 2-D PCs involving anisotropic materials. We will investigate a structure involving liquid crystals (LCs). The LC is a uniaxial material and its optical properties depend on the direction of the propagation and polarization state of lightwaves relative to the orientation of the LC director. When the director is perpendicular to the PC extension, the E and H waves are decoupled. Here, the E (H) wave refers to the mode with the electric (magnetic) field parallel to the PC extension direction. Therefore, the E wave is like that in an isotropic PC. We will study the variation of band gaps for the H wave for different radii of the LC cylinders and different orientations of the LC director in square lattice, and investigate the symmetry properties of the PBG with respect to the orientation of the LC director.

2. The FDFD formulation

The FDFD formulation starts from the six component equations of the two Maxwell curl equations for a general anisotropic dielectric medium. The z-axis is defined as along the extension direction of the 2-D PC. Consider the situation that one principal axis of the anisotropic medium is along the z direction such that the permittivity tensor elements ²xz, ²yz, ²zx, and ²zyare zeros. We use the Yee-mesh-based finite-difference

scheme [5]. For the case of H wave, since the different components of the electric field intensity vector E and the electric flux density vector D are assigned at different grid points in the Yee mesh, the tensor constitutive relation between E and D must be carefully treated at a given point. We employ a scheme of making spatial averaging of the suitable component of D, as given by the following expressions [6]:

Ex(i + 0.5, j) =²yy Λ Dx(i + 0.5, j) − ²xy 4Λ{Dy(i, j − 0.5) + Dy(i + 1, j − 0.5) + Dy(i, j + 0.5) + Dy(i + 1, j + 0.5)} (1a)

Appendix II

Ey(i, j + 0.5) = − ²yx

4Λ{Dx(i + 0.5, j) + Dx(i + 0.5, j + 1) + Dx(i − 0.5, j) + Dx(i − 0.5, j + 1)} +²xx

Λ Dy(i, j + 0.5) (1b)

where Λ = ²0(²xx²yy− ²xy²yx) and ²0is the free-space permittivity. Finally, an eigenvalue equation in terms

of Hz in the following form can be obtained:

k20Hz= − {UxB22Vx− UxB21Vy+ VyA11Uy− UyA12Vx}Hz (2)

where Ux, Uy, Vx, Vy, B22, B21, A11, and A12 are all matrices. For the E wave, the eigenvalue equation is

obtained as

k2

0Ez= − ²−1zz{VxUx+ VyUy}Ez. (3)

In the analysis of PCs, the periodic boundary conditions (PBCs) need to be considered with the unit cell. For a square unit cell with width a, the PBCs for the field ψ to be solved can be expressed as

ψ(x, y + a) =e−jkya_{ψ(x, y) (4a)}

ψ(x + a, y) =e−jkxa_{ψ(x, y). (4b)}

3. Numerical results

As a numerical example, we consider a 2D PC composed of circular cylinder LC-filled regions with silicon

Figure 1: (a) The cross-section of the square-lattice 2D PC considered. (b) The entire first Brillouin zone. as the background medium. The refractive index of silicon is taken to be n = 3.4. Figure 1(a) shows the cross-section of this 2D PC, where the cylinder rod of radius r = 0.5a is the region with the nematic LC, where a is the period of the lattice. The nematic LC is a homogenous uniaxial LC having the ordinary and extraordinary refractive indices no = 1.5292 and ne = 1.7072, respectively. The permittivity tensor

elements in the LC region are expressed as

²xx=n2o+ (n2e− n2o) sin2θ cos2α (5a) ²xy=²yx= (n2e− n2o) sin2θ sin α cos α (5b) ²yy=n2o+ (n2e− n2o) cos2θ sin2α (5c)

where θ is the angle between the LC director and the z axis, and α is the angle between the projection of the LC director on the x-y plane. We assume θ = 90◦ _{in this example. The first Brillouin zone (BZ)}

of the square lattice is shown in Fig. 1(b). With the existence of anisotropic materials, the symmetry properties of the band structure become more complicated. For the present PC, we need to consider half of the first BZ. Figure 2 presents the band structure for the PC with α = 45◦ _{for the H wave. This result}

has been checked with the calculation conducted by the finite element method [7] and very good agreement is achieved. A bandgap exists between the fourth and fifth bands. The calculated gap-to-midgap ratio versus the cylinder radius for α = 45◦ _{is shown in Fig. 3(a). It is seen that when r = 0.5a, the maximum}

Figure 2: The band structure for α = 45◦

radius (lattice constant)

( )a ( )b a (degree) 0.4 0.42 0.44 0.46 0.48 0.5

a=45

Figure 3: (a) Gap-to-midgap ratio as a function of the cylinder radius. (b) Gap-to-midgap ratio as a function of α.

gap size occurs. The calculated gap-to-midgap ratio as a function of α is plotted in Fig. 3(b), which shows symmetry with respect to 45◦_.

4. Conclusion

We have developed an Yee-mesh-based FDFD method for calculating the band structures of 2-D PCs containing anisotropic materials. As a numerical example, we have analyzed a square-lattice PC composed of circular cylinders of nematic LC material with silicon as background. The LC director is assumed to be perpendicular to the PC extension direction so that the E and H waves are decoupled. The tunability of the band-gap size has been studied and symmetry properties of the gap size versus the orientation of the LC director observed. The FDFD calculation has been found to agree with the finite element analysis [7]. 5. References

[1]E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. 58, 2059–2062 (1987).

[2]S. John, “Strong localization of photons in certain disordered dielectric superlattices,” Phys. Rev. Lett. 58, 2486–2489 (1987).

[3]I. H. H. Zabel and D. Stround, “Photonic band structures of optically anisotropic periodic arrays,” Phys. Rev. B 48, 5004–5012 (1993).

[4]Z.-Y. Li, J. Wang, and B.-Y. Gu, “Creation of partial bandgaps in anisotropic photonic-band-gap structures,” Phys. Rev. B 12, 1397–1408 (2004).

[5]C. P. Yu and H. C. Chang, “Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals,” Opt. Express 58, 3721–3729 (1998).

[6]A. P. Zhao, J. Juntunen and A. V. Raisanen, “An efficient FDTD algorithm for the analysis of microstrip patch antennas printed on a general anisotropic dielectric substrate,” IEEE Trans. Microwave Theory Tech. 47, 1142–1146 (1999). [7] S. M. Hsu, M. M. Chen, and H. C. Chang, “Band-structure analysis of 2D non-diagonal anisotropic photonic crystals by the finite element method,” submitted to IPNRA 2007.

Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method

Po-jui Chiang,1Chin-ping Yu,1and Hung-chun Chang1,2,3,*

1_{Graduate Institute of Electro-Optical Engineering, National Taiwan University, Taipei, Taiwan 106-17, Republic of China}

2

Graduate Institute of Communication Engineering, National Taiwan University, Taipei, Taiwan 106-17, Republic of China

3

Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan 106-17, Republic of China

共Received 18 April 2006; revised manuscript received 21 November 2006; published 20 February 2007兲

An analysis method based on a multidomain pseudospectral method is proposed for calculating the band diagrams of two-dimensional photonic crystals and is shown to possess excellent numerical convergence behavior and accuracy. The proposed scheme utilizes the multidomain Chebyshev collocation method. By applying Chebyshev-Lagrange interpolating polynomials to the approximation of spatial derivatives at collo-cation points, the Helmholtz equation is converted into a matrix eigenvalue equation which is then solved for the eigenfrequencies by the shift inverse power method. Suitable multidomain division of the computational domain is performed to deal with general curved interfaces of the permittivity profile, and field continuity conditions are carefully imposed across the dielectric interfaces. The proposed method shows uniformly ex-cellent convergence characteristics for both the transverse-electric and transverse-magnetic waves in the analy-sis of different structures. The analyanaly-sis of a mini band gap is also shown to demonstrate the extremely high accuracy of the proposed method.

DOI:10.1103/PhysRevE.75.026703 PACS number共s兲: 02.70.Hm, 03.50.De

I. INTRODUCTION

Band structures are essential characteristics of photonic crystals 共PCs兲, from which possible photonic band gaps 共PBGs兲 can be identified 关1–3兴. For frequencies within the

PBGs, wave propagation is forbidden and many photonic devices have been proposed and designed based on this phe-nomenon. In particular, two-dimensional 共2D兲 PCs com-posed of either dielectric rods or air columns have been widely employed in many applications such as waveguiding, resonant cavity formation, and wavelength filtering. In this paper, we propose an analysis scheme with excellent numeri-cal convergence behavior and accuracy for numeri-calculating the band structures of 2D PCs. The currently most used numeri-cal methods for such numeri-calculations have been the plane-wave expansion 共PWE兲 method 关3–6兴 and the finite-difference

time-domain共FDTD兲 method 关7,8兴. The finite-difference

ei-genvalue problem formulation has also been employed by Yang关9兴 and Shen et al. 关10兴, and more recently by Yu and

Chang关11兴 based on the Yee mesh as often employed in the

FDTD method 关12兴. The Yee-mesh-based formulation was

named the finite-difference frequency-domain 共FDFD兲 method. Yu and Chang关11兴 used the FDFD method to

ana-lyze the band structures of 2D PCs with either square or triangular lattices and adopted a fourth-order accurate com-pact finite-difference scheme关13兴 to increase numerical

effi-ciency and accuracy. Although the FDFD method offers re-sults with accuracy comparable to those obtained using the MIT photonic-bands共MPB兲 package 关14兴 based on the PWE

method, the numerical convergent speed was found not to be uniformly fast among different bands in the two methods.

The numerical formulation proposed in this paper is based on the multidomain pseudospectral method using Chebyshev polynomials. The pseudospectral method has recently

at-tracted raised attention as an alternative treatment for com-putational electromagnetics because of its high-order accu-racy and fast convergence behavior over traditional techniques while retaining formulation simplicity. It has a long history of being applied to fluid dynamics关15兴 and has

recently been extended to the analysis of electromagnetics both in the time关16–18兴 and in the frequency domain 关19兴.

However, while the theory of the pseudospectral method has been well elaborated, the application in the frequency do-main has not received much focus compared with that in the time domain in the electromagnetics community. In关19兴, the

pseudospectral frequency-domain method was proposed and applied to solve the nonhomogeneous 共nonzero-source兲 Helmholtz equation in a simple two-subdomain problem with rectangular structure shape. Our proposed scheme in this paper utilizes the multidomain Chebyshev collocation method supported by the curvilinear mapping technique关20兴

to facilitate and ameliorate the simulation of 2D PCs of ar-bitrary permittivity profile. The formulation is derived in the form of an eigenvalue problem so that we can readily obtain the eigenmodes by available mathematical tools. Here, we adopt the shift inverse power method共SIPM兲 for its particu-larly fast convergence characteristic over other conventional methods applying matrix inversion. Both the multidomain pseudospectral algorithm and the SIPM furnish our pseu-dospectral mode solver共PSMS兲 as a quite powerful and flex-ible method. To obtain high-accuracy full-vectorial modal solutions for dielectric structures, proper satisfaction of di-electric interface conditions is essential, whether it is based on the finite-difference method 关21,22兴, the finite-element

method 关23兴, or others. Such proper treatment of interface

conditions will be carefully considered in our formulation. The rest of this paper is outlined as follows. The physical problem involving the Helmholtz equations is described in Sec. II along with the required Dirichlet and Neumann type boundary conditions across the dielectric interfaces. The for-mulation of the PSMS is presented in Sec. III. Numerical examples including 2D PCs with either a square or a

trian-*Electronic address: [email protected]

PHYSICAL REVIEW E 75, 026703共2007兲

1539-3755/2007/75共2兲/026703共14兲 026703-1 ©2007 The American Physical Society

NEAR-FIELD AND FAR-FIELD BEHAVIOR OF A METALLIC

NANOSCALE SPHERE AT OPTICAL FREQUENCIES BASED ON THE

CLASSICAL DRUDE MODEL

Bang-Yan Lin

1

,

Chun-Hao Teng

2

, and Hung-chun Chang

1

1

_{Dept. of Electrical Eng., Graduate Institute of Electro-Optical Eng., and Graduate Institute of}

Comm. Eng., National Taiwan Univ., Taipei, Taiwan 106-17, R.O.C.

[email protected]

2

_{Math. Dept., National Cheng Kung Univ., Tainan, Taiwan 701, R.O.C.}

_{[email protected]}

Abstract: Optical properties of a metallic sphere based on the classical Drude model are studied theoretically. The

comparison between the metallic sphere and the perfectly conducting one shows conspicuous distinction in the electric field amplitude on the surface of the sphere and in the radar cross section behavior.

INTRODUCTION

Recently, the field of nano-optics has been rapidly developing due to mature fabrication

technologies. One relevant important subject is the phenomenon of the surface plasmon that was

induced by the interaction of metal with electromagnetic radiation. The phenomenon is

especially evident in some novel metals, such as silver or gold. In the scattering measurement

problem, one of the usual geometric structures is the particle that can be approximated by a

sphere. Fortunately, there has existed a well-known analytic solution for the spherical-particle

scattering, i.e., the Mie solution [1], which not only can be applied to a perfectly conducting

sphere but also to a metallic one with a complex value of the permittivity. The theory provides an

initial prediction of the behavior in the near and far regions of the sphere. Traditionally, the use

of the perfectly conducting sphere in the scattering problem can provide good approximation in

the microwave region. But when the frequency moves to the optical regime, such as the visible

region, and the particle shrinks to the nanometer-scale, one should care the appropriateness of

using the perfectly conducting sphere. In this paper, to conform the real situation, we consider

the classical Drude model for a metallic sphere. We discuss the near field and far field behavior,

respectively. The study of the near field is aimed at whether the size of the particle will lead to a

significant difference in the amplitude of electric field on the surface of the sphere. As for the far

field, the radar cross section (RCS) will be reconsidered under the new situation. Besides using

the analytic solution, we also employ a high-accuracy 3-D time-domain numerical scheme to

verify the RCS results.

NEAR-FIELD BEHAVIOR

Fig. 1 shows a plane wave with the x-direction polarization propagating along the z direction and

incident on a metallic sphere, where the dielectric constants of the background space and the

Fig. 1. Uniform plane wave incident on a metallic sphere.

sphere are

_ε

( )I

_and

_ε

( )II

_{, respectively, and the radius of the sphere is a . We assume}

( )

1.0

I

ε

=

Appendix IV

and the dielectric constant of the sphere is dominated by the classical Drude model with its

dispersion characteristic expressed as

( )II

_{( )}

2

_/(

2

_{/ ).}

p

i

ε

ω

=

ε

∞

−

ω ω

+

ω τ

Here

ε

∞

is the

high-frequency limit of the dielectric constant,

ω

_p

is the plasma frequency, and

τ

is the

relaxation time. Fig. 2 shows the experimental data of the dielectric constant of silver referring to

[2], and the curve-fitting results based on two different sets of parameters for the classical Drude

model, the first one being

_{4.0, 133.6848 10 Hz,}

14

_{and 0.9019 10}

14

_sec

p

ε

ω

τ

−

∞

=

=

×

=

×

from [3]

and the second one being

_{3.7, 136.73 10 Hz,}

14

_{and 3.6563 10 sec}

14 p

ε

ω

τ

−

∞

=

=

×

=

×

from [4],

denoted as classical Drude models 1 and 2, respectively.

(a) (b)

Fig. 2. (a) The real part and (b) the imaginary part of the dielectric constant of silver.

Fig. 3 shows the maximum amplitude of the electric field on the surface of the sphere in the

wavelength range

λ

=

0.25 μm ~ 1 μm

, where “Experimental Data” mean calculated results

using the corresponding

_ε

( )II

_{shown in Fig. 2. Fig. 3 (a) and (b) is for}

_a

₌

_{0.25 μm}

_{and Fig. 3}

(c) and (d) for

a

=

0.025 μm

, with Fig. (b) and (d) being the expanded versions for

0.25 μm ~ 0.45 μm

λ

=

We can see that the values for the silver sphere are larger than those for

the perfectly conducting one for most of the wavelengths. In particular, for

λ

=

0.3 ~ 0.4 μm

,

the differences are more significant. In both cases, the amplitude reaches the highest

for

λ

:

0.36 μm

, partly due to the minimum in the imaginary part shown in Fig. 2(a). The highest

values represented by the experimental data are about 7 and 19, respectively. The results show

that the size of the sphere will affect the amplitude of the electric field and that for the metallic

sphere is larger than that of the perfectly conducting one.

(a) (b) (c) (d)

Fig. 3. The maximum amplitude of the electric field on the surface of the sphere. The radius of the sphere is

0.25 μm

in (a) and (b). The radius of the sphere is

0.025 μm

in (c) and (d).

FAR-FIELD BEHAVIOR

The measurement of electromagnetic scattering in the far-field zone is one of the important

approaches for investigating the properties of the target, which is well known to be described by

the RCS quantitatively. In our study, we use the classical Drude model 1 because it fits better the

experimental data than the type 2. The definition of RCS can be written as

2 2 2 s i 3-D

( , , ) lim[4

_r

r

E

/ E ]

σ

λ θ φ

π

→∞

=

,

where

_{E is the incident field and}

i

_{E is the scattering}

s

field. Here we discuss two kinds of RCS, bistatic RCS and monostatic RCS, with calculated

results shown in Fig. 4 (a) and (b), respectively. The parameters for Fig. 4(a) are

λ

=

0.364 μm

,

[0, ]

θ

=

π

, and

φ π

= , and those for Fig. 4(b) are /

a

λ

=

0 ~ 1.5

,

θ π

=

, and

φ π

= . Both results

are normalized to

_π

_a

2

_{and a used in the calculation is}

_{0.25 μm . Besides comparing the exact}

results of the silver sphere with those of the perfectly conducting one, we also conduct

time-domain numerical simulations using a high-accuracy scheme similar to that of [5] along

with the technique of near-to-far field transformation [6]. We can see that the differences

between the metallic sphere and the perfectly conducting one are obvious in both cases except in

the Rayleigh region (

a

<

0.1

λ

) for the monostatic RCS. The numerical results for the bistatic

RCS agree very well with the exact results. But for the monostatic RCS, there exists deviation at

high frequencies. This is attributed to lack of enough resolution. Using denser meshes in the

numerical scheme can improve the results.

(a) (b)

Fig. 4. (a) Normalized bistatic RCS. (b) Normalized monstatic RCS.

CONCLUSION

At optical frequencies, we have discovered that the size of silver sphere will make a large

distinction on the amplitude of the electric field in the near field, as predicted by the Mie theory.

In the far field, the RCS properties of the metallic sphere are unlike those of the perfectly

conducting one. These theoretical results could provide a useful reference for related

experiments involving nanoscale metallic spheres.

REFERENCES

[1] Born, M., and Wolf E. “Principle of Optics”, Cambridge University Press, U.K., 1999.

[2] Johnson, P. B., and Christy R. W., “Optical Constants of Noble Metals,” Phys. Rev. B., vol. 6,

no. 12, Dec. 1972, pp. 4370–4379.

[3] Etchegoin, P. G., and Ru L. E. C., “Multipolar emission in the vicinity of metallic

nanostructures,” J. Phys.: Condens. Matter, vol. 18, 2006, pp. 1175–1188.

[4] Saj, W. M., “FDTD simulations of 2D plasmon waveguide on silver nanorods in hexagonal

lattice,” Opt. Express, vol. 13, 2005, pp. 4818–4827.

[5] Hesthaven, J. S., and Gottlieb D., “A stable penalty method for the compressible

Navier-Stokes equations. III. Multidimensional domain decomposition schemes,” SIAM J. Sci.

Comput., vol. 20, 1998, pp. 62–93.

[6] Taflove, A., and Hagness A. C., “Computational Electrodynamics: The Finite-Difference

Air-Core Waveguides for Terahertz Signals

Chin-ping Yu and Hung-chun Chang*

Graduate Institute of Electro-Optical Engineering, National Taiwan University, Taipei, Taiwan 106-17, R.O.C.

*also with Department of Electrical Engineering and Graduate Institute of Communication Engineering, National Taiwan University E-mail: [email protected]

Abstract: The finite-difference frequency-domain method is adopted for the analysis of air-core

THz waveguides formed by Teflon rods. The guiding mechanism is found to be based on the ARROW model. The calculated propagation characteristics demonstrate that well confined THz guided field can be obtained.

©2007 Optical Society of America

OCIS codes: (060.2310) Fiber optics; (230.7370) Waveguides; (260.2110) Electromagnetic theory; (999.9999) THz optics.

1. Introduction

In most terahertz (THz) systems, the propagation of far infrared electromagnetic waves relies on traditional metal or dielectric waveguides which suffer from high conductivity losses or dielectric losses and consequently have limitation for practical applications. Thus, it is a great issue to develop effective and low-loss waveguiding structures for THz signals for a compact, reliable, and flexible THz system. One possible way is to use metal-based Sommerfeld wires to confine THz signals with low fraction of power within lossy dielectric materials to reduce the attenuation caused by materials [1]. However, the transmission of THz signals on Sommerfeld wires suffers from a large field extension into air and large radiation loss at waveguide bending. Another approach utilizes total internal reflection (TIR) to guide THz signals in small air region formed by dielectric materials [2]. Good confinement of THz signals can be observed in the small air gap but with unignorable amount of field in the thick dielectric layers, which might result in high transmission loss in the materials in THz frequency range.

In this paper we consider air-core waveguides formed by polytetrafluorethylene (Teflon) rods for THz waveguiding. Figs. 1(a) and 1(b) show two examples of such waveguides with the air core surrounded by variant numbers of Teflon rods. The propagation characteristics of the waveguides with various rod-sizes and distributions are calculated by applying a full-vectorial finite-difference frequency-domain (FDFD) method [3]. The performance of these THz waveguides is discussed based on the calculated propagation characteristics.

i+1/2 j+1/2 Ex(i-1/2, j ) Ey( i , j-1/2) Ez(i , j) ε1 ε2 H_x( i , j-1/2) H_z(i-1/2, j-1/2) Hy(i-1/2, j ) Hz(i +1/2, j+1/2) i-1/2 j-1/2 ε_z Hz(i +1/2, j-1/2) Hz(i -1/2, j+1/2) Ex(i+1/2, j ) Hy(i+1/2, j ) Ey( i , j+1/2) Hx( i , j+1/2)

Fig. 1. (a) The cross-section of an core waveguide with 8 Teflon rods. (b) The cross-section of an air-core waveguide with 12 Teflon rods. (c) Yee’s mesh for the FDFD method.

2. The numerical method

The FDFD method which is an efficient and accurate numerical model for optical waveguide analysis [3] is utilized to analyze air-core waveguides in Figs. 1(a) and 1(b). Consider the time-harmonic Maxwell's curl equations

H j E =− ωμ₀ ×

∇ and ∇×H = jωε0εrE where ω is the angular frequency, μ0 and ε0 are the permittivity and the

permeability of free space, respectively, and εr is the relative permittivity of the dielectric medium. By discretizing

the transverse (x-y) plane using the Yee’s mesh as shown in Fig. 1(c) and applying the central difference scheme for the differential operators, the curl equations can be converted into the matrix form [3]

(1) ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − − = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − − − = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − z y x x y x y z y x z y x z y x x y x y z y x H H H 0 V V V 0 I V I 0 E E E ε ε ε E E E U U U I U I H H H jβ jβ j j j j 0 0 0 0 0 0 0 0 0 0 0 ωε β β ωμ

where β is the modal propagation constant, H and E are the vectors composed of the electric and magnetic field components, respectively, at the grid points, I is a square identity matrix, and εi's (i = x, y, z) are diagonal matrices

representing the relative permittivities at the corresponding grid points. Ux, Uy, Vx, and Vy are square matrices

determined by the central difference scheme and the boundary conditions. After some mathematical work, an eigenvalue matrix equation in terms of the transverse magnetic fields can be obtained as

2 _. β ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ x xx xy x x y yx yy y y H Q Q H H Q H Q Q H H (2)

Solving the eigenvalue problem of (2) by the shift inverse power method, the propagation constants and the field distributions of guided modes on the air-core waveguides can be obtained.

3. Results and discussion

In the following computation, the perfectly matched layers (PMLs) are applied as the boundary conditions around the computational window and the material loss is taken into account. The real part of the refractive index of Teflon is set to be 1.443 independent of frequency, and the value of its imaginary part is obtained by fitting the experimental results [4] with a second order polynomial as a function of frequency. We first consider an air-core waveguide with the air-core surrounded by 8 Teflon rods as shown in Fig. 1(a). The core size R is set to be 4.0 mm and the diameter of Teflon rods is d = 0.3R = 1.2 mm. By applying the FDFD method within the frequency range from 750 GHz to 1100 GHz, the modal refractive indices and the corresponding propagation losses of guided modes on such waveguide are shown in Figs. 2(a) and 2(b), respectively. The propagation loss is obtained from the modal propagation constant by the relation: Loss (dB/m) = 20·103_{·Im(β)/ln(10). In Fig. 2(a) it can be observed that over the}

entire frequency range, there always exist guided modes except at some particular frequencies. Fig. 2(c) illustrates the guided magnetic field profiles at f = 850 and 980 GHz. At f = 850 GHz, one can see that most of the field is well confined in the air core with relatively small field in the Teflon rods, thus reducing the possible attenuation loss caused by materials. As the frequency approaches the discontinuities in the modal dispersion curve, relatively large field appears in the Teflon rods and penetrates into the outer air region, as shown in Fig. 2(c) for f = 980 GHz, resulting in relatively large propagation losses at such frequencies shown in Fig. 2(b). The guiding mechanism of the air-core waveguide is just like the antiresonant reflecting optical waveguide (ARROW) model discussed in [5]. The discontinuities along the modal dispersion curve correspond to resonant frequencies of the Teflon rods at which the mode field leaks into the cladding region as in Fig. 2(c). At other frequencies, the field can be reflected by the Teflon rods to support centrally localized guided modes with relatively low propagation losses.

7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 0.9988 0.9990 0.9992 0.9994 0.9996 0.9998 M odal index Frequency (Hz) x1011 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 L oss (d B /m ) Frequency (Hz) x1011

Fig. 2. (a) The modal index and (b) the propagation loss of guided modes on the air-core waveguide of Fig. 1(a). (c) Field distributions at f = 850 and 980 GHz.

If we keep the size of each Teflon rod the same and increase the number of rods surrounding the air core, similar propagation characteristics of guided modes can be obtained. Fig. 3(a) shows the comparison of modal dispersion curves for 8 and 12 Teflon rods. Slightly smaller modal indices are observed for 12 Teflon rods, and the

positions of discontinuities are almost the same for both arrangements. Besides, the propagation-loss results given in Fig. 3(b) also show similar distributions for 8 and 12 rods with the peaks appearing at similar locations.

7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 8 rods 12 rods L oss ( dB /m ) Frequency (Hz) x1011 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 0.9988 0.9990 0.9992 0.9994 0.9996 0.9998 8 rods 12 rods M odal index Frequency (Hz) x1011

Fig. 3. (a) The modal index and (b) the propagation loss of guided modes on air-core waveguides with variant numbers of rods.

Now we keep the size of the air core the same and change the diameter of the Teflon rods. Figs. 4(a) and 4(b) show the modal indices and propagation losses, respectively, for 8 teflon rods with d = 1.2 mm and 1.6 mm. One can see that the discontinuities in modal dispersion curves and the positions of the propagation-loss peaks have obvious differences between the two sizes. This is because the resonant frequencies highly dependent on the size of the surrounding rods and so do the propagation characteristics. If the sizes of the Teflon rods are kept the same, the number of rods will not affect the resonant frequencies, making the air-core waveguides possess similar propagation characteristics as shown in Fig. 3. Practically, these rods need an inner or outer shell to be connected to realize a workable structure, which might influence the propagation properties and is worth further discussions.

7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 R = 4.0 mm d = 1.2 mm R = 4.0 mm d = 1.6 mm L oss (d B /m ) Frequency (Hz) x1011 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 0.9988 0.9990 0.9992 0.9994 0.9996 0.9998 R = 4.0 mm d = 1.2 mm R = 4.0 mm d = 1.6 mm M od al in de x Frequency (Hz) x1011

Fig. 4. (a) The modal index and (b) the propagation loss of guided modes on 8-rod air-core waveguides with variant size of rods.

4. Conclusions

We have investigated modal properties of guided modes on air-core THz waveguides with the air core surrounded by Teflon rods. The guided modes can be found over almost the frequency range concerned and the ARROW model is successfully adopted to explain the discontinuities of the modal index curves where the propagation losses become relatively huge. It is also observed that the propagation characteristics of the air-core waveguides are highly dependent on the size of surrounding rods and the number of rod has little influence on the performance of these waveguides.

5. References

[1] K. Wang and M. Mittleman, “Metal wires for terahertz wave guiding,” Nature 432, 376–379 (2004).

[2] M. Nagel, A. Marchewka, and H. Kurz, “Low-index discontinuity terahertz waveguides,” Opt. Express 14, 9944–9954 (2006).

[3] C. P. Yu and H. C. Chang, “Yee-mesh-based finite difference eigenmode solver with PML absorbing boundary conditions for optical waveguides and photonic crystal fibers,” Opt. Express 12, 6165–6177 (2004).

[4] J. R. Birch, J. D. Dromey, and J. Lesurf, “The optical constants of some common low-loss polymers between 4 and 40 cm-1_{,” Infrared Phys.} 21, 225–228 (1981).

[5] T. P. White, R. C. McPhedran, C. M. Sterke, N. M. Litchinitser, and B. J. Eggleton, “Resonance and scattering in microstructured optical fibers,” Opt. Lett. 27, 1977–1979 (2002).

JOURNAL OF LIGHTWAVE TECHNOLOGY, VOL. 24, NO. 11, NOVEMBER 2006 4411

Robust Calculation of Chromatic

Dispersion Coefficients of Optical Fibers From

Numerically Determined Effective Indices Using

Chebyshev–Lagrange Interpolation Polynomials

Po-Jui Chiang, Chin-Ping Yu, and Hung-Chun Chang, Senior Member, IEEE, Member, OSA

Abstract—Numerical calculation of chromatic dispersion

co-efficients of optical fibers is conducted using a procedure involving Chebyshev–Lagrange interpolation polynomials. Only numerically determined effective indices at several wavelengths are needed for obtaining the dispersion curve, and no direct nu-merical differentiation of the effective refractive index is involved. A silica-filled metallic rectangular waveguide having analytical solutions for the effective refractive index and the chromatic dispersion is used as an example for confirming the accuracy and efficiency of the proposed method. The method is then also applied to the analysis of holey fibers.

Index Terms—Dispersion, fiber properties, optical waveguides,

photonic-crystal fibers (PCFs), spectral method.

I. INTRODUCTION

T

HE CHROMATIC dispersion coefficient of an optical fiber in units of picoseconds per nanometer kilometer is a key quantity in various analysis and design issues in fiber transmission systems, including recent intensive investigation of photonic-crystal fibers (PCFs) or holey fibers [1]–[3]. It is defined as D= −(λ/c)d2n_eff/dλ2, where n_eff is the effective refractive index of the optical fiber mode, λ is the wavelength, and c is the speed of light in a vacuum. For a given fiber or waveguide structure, neffcan be solved at different wavelengths using various numerical methods. Then, in principle, one sim-ple way to obtain D at a particular wavelength is to perform direct numerical differentiation using the central difference scheme based on three n_eff values at three nearby wavelengths. Therefore, to obtain D values at N wavelengths, one needs to calculate n_eff at3N wavelengths. From the basic theory of

Manuscript received May 27, 2006; revised August 14, 2006. This work was supported in part by the National Science Council of the Republic of China under Grant NSC94-2215-E-002-022.

P.-J. Chiang and C.-P. Yu are with the Graduate Institute of Electro-Optical Engineering, National Taiwan University, Taipei 106-17, Taiwan, R.O.C. (e-mail: [email protected]; [email protected]).

H.-C. Chang is with the Department of Electrical Engineering, Gradu-ate Institute of Electro-Optical Engineering, and the GraduGradu-ate Institute of Communication Engineering, National Taiwan University, Taipei 106-17, Taiwan, R.O.C. (e-mail: [email protected]).

Color versions of Figs. 1, 3, 4, and 6 are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/JLT.2006.883646

the finite difference method, the separation between adjacent wavelengths should be small enough to maintain the required accuracy but should not be too small such that the error would explode [4]. Thus, some suitable ∆λ has to be determined. Furthermore, one might worry about the effect of the numerical accuracy of the calculated neff on the degree of accuracy of its second derivative. Kuhlmey et al. [5] have reported in their analysis of holey fibers that they generally considered ∼1000 points per unit λ/Λ, with Λ being the pitch of the air-hole lattice of the air-holey fiber, to obtain satisfactory results of dispersion curves.

To ensure accurate determination of D, over the past two decades, some more elaborate approaches have been proposed including, for example, that based on the matrix perturba-tion method [6] and those through formulating two associated problems for evaluating the first and second derivatives of the propagation constant [7]–[10]. In this paper, we propose a simple approach utilizing Chebyshev–Lagrange interpola-tion polynomials as in the Chebyshev spectral method or the Chebyshev collocation method (CCM) [11] and demonstrate its efficiency and robustness in calculating D. For instance, over the wavelength range from 0.6 to 1.5 µm, calculation of neff at only 13 wavelength points would guarantee accurate prediction of the D curves. Most importantly, no numerical differentiations of n_eff are required, and thus, computational burden can be greatly reduced.

The rest of this paper is organized as follows: The CCM is described in Section II. The proposed procedure for chro-matic dispersion coefficient calculation is given in Section III, along with the analysis of a silica-filled metallic rectangular waveguide. The application to holey fibers and comparison with published results are presented in Section IV. Conclusions are drawn in Section V.

II. CCM

We first explain the approximation of a function f(x) in the domain−1 ≤ x ≤ 1 under the CCM [11]. The Chebyshev polynomial TN(x) of degree N is defined as TN(x) =

cos(N cos−1_{x), where |x| ≤ 1, and the collocation points to}

be used are given by the Chebyshev–Gauss–Lobatto points defined as the roots of the polynomial (1 − x2)T_N , where

0733-8724/$20.00 © 2006 IEEE