Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions

doi:10.1364/OE.14.011631

« Show journal navigation

Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions

Chia-Chien Huang »View Author Affiliations

Optics Express, Vol. 14, Issue 24, pp. 11631-11652 (2006)
http://dx.doi.org/10.1364/OE.14.011631

View Full Text Article

Acrobat PDF (760 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Formulations of planar waveguides
3. Numerical scheme
4. Numerical results and discussion
5. Conclusion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (26)
Cited By
Figures (9)
Metrics

Abstract

The pseudospectral method, proposed in our previous work, has not yet been constructed for optical waveguides with leaky modes or anisotropic materials. Our present study focuses on antiresonant reflecting optical waveguides (ARROWS) made by anisotropic materials. In contrast to the fields in the outermost subdomain expanded by Laguerre-Gaussian functions for guided mode problems, the fields in the high-index outermost subdomain are expanded by the Chebyshev polynomials with Mur’s absorbing boundary condition (ABC). Accordingly, the traveling waves can leak freely out of the computational window, and the desirable properties of the pseudospectral scheme, i.e., provision of fast and accurate solutions, can be preserved. A number of numerical examples tested by the present approach are shown to be in good agreement with exact data and published results achieved by other numerical methods.

1. Introduction

Building optical integrated circuits for practical use in an optical communication system requires a variety of photonic components such as optical interconnectors, modulators, switches, filters, splitters, polarizers, and couplers. The diverse optical waveguides by using distinct materials and refractive index profiles, according to specifications, produce functional optical integrated circuits. To precisely compute propagation characteristics (i.e., dispersion, attenuation, mode distribution, and birefringence) of general optical waveguides with both guided and leaky modes, we may precisely tailor the expected optical devices. A low-loss waveguide antiresonant reflecting optical waveguide (ARROW), which utilizes two thin interference cladding films made on semiconductor substrates [1

M. A. Duguay, Y. Kokubun, T. L. Koch, and L. Pfeiffer, “Antiresonant reflecting optical waveguides in SiO₂-Si multilayer structures,” Appl. Phys. Lett. 49, 13–15 (1986). [CrossRef]

], has been proposed to construct a wide range of functional devices [2–5

T. Baba, Y. Kokubun, T. Sakaki, and K. Iga, “Loss reduction of an ARROW waveguide in shorter wavelength and its stack configuration,” J. Lightwave Technol. 6, 1440–1445 (1988). [CrossRef]

]. Unlike total internal reflection mechanisms for guided modes, ARROW structure supports leaky waves through high reflectivity corresponding to the antiresonant condition of a Fabry-Perot resonator. This kind of waveguide has several advantages [6

T. Baba and Y. Kokubun, “Dispersion and radiation loss characteristics of antiresonant reflecting optical waveguides-numerical results and analytical expressions,” IEEE J. Quantum Electron. 28, 1689–1700 (1992). [CrossRef]

W. P. Huang, R. M. Shubair, A. Nathan, and Y. L. Chow, “The modal characteristics of ARROW structures,” J. Lightwave Technol. 10, 1015–1022 (1992). [CrossRef]

] such as low loss, relatively large spot size for effectively coupling into optical fibers, high fabrication tolerance for thicknesses of interference claddings, high loss discrimination between the fundamental and higher-order modes, high polarization sensitivity, and polarization insensitivity [8

T. Baba and Y. Kokubun, “New polarization-insensitive antiresonant reflecting optical waveguide (ARROW-B),” IEEE Photon. Technol. Lett. 1, 232–234 (1989). [CrossRef]

In the viewpoint of computation, ARROW-type devices may be regarded as multilayer waveguide structures, and thus various analytical or numerical approaches including transverse resonance method (TRM) [7

W. P. Huang, R. M. Shubair, A. Nathan, and Y. L. Chow, “The modal characteristics of ARROW structures,” J. Lightwave Technol. 10, 1015–1022 (1992). [CrossRef]

W. Jiang, J. Chrostowski, and M. Fontaine, “Analysis of ARROW waveguides,” Opt. Commun. 72, 180–186 (1989). [CrossRef]

] and transfer matrix method (TMM) [10–12

J. Chilwell and I. Hodgkinson, “Thin-film field-transfer matrix theory for planar multilayer waveguides and reflection from prism-loaded waveguides,” J. Opt. Soc. Am. A 1, 742–753 (1984). [CrossRef]

] have been developed to find the dispersion equations. To calculate propagation characteristics, except some traditional zero-search algorithms like Newton’s and scan methods, which cannot find all of the guided and leaky modes simultaneously, a more efficient and reliable technique is proposed for extracting the poles of dispersion equations in complex plane by the argument principle method (APM) [13

E. Anemogiannis and E. N. Glytsis, “Multilayer waveguides: efficient numerical analysis of general structures,” J. Lightwave Technol. 10, 1344–1351 (1992). [CrossRef]

,14

C. K. Chen, P. Berini, D. Feng, S. Tanev, and V. P. Tzolov, “Efficient and accurate numerical analysis of multilayer planar waveguides in lossy anisotropic media,” Opt. Express 7, 260–272 (2000). [CrossRef] [PubMed]

]. However, the TRM and TMM are only accurate enough to be used with planar waveguides that have step-index profiles rather than graded-index profiles, because inhomogeneous materials are often approximated by a stack of dielectric layers with constant refractive indices. Moreover, the two schemes mentioned above are also unsuitable for the 2-D cross-section structures like rib ARROW providing also lateral confinement on mode profiles. Accordingly, it is necessary and desirable to pay close attention to the rigorous numerical solution methods. The schemes used most often such as the finite difference modal analysis [15

W. P. Huang, C. L. Xu, W. Lui, and K. Yokoyama, “The perfectly matched layer boundary conditions for modal analysis of optical waveguides: leaky mode calculations,” IEEE Photon. Technol. Lett. 8, 652–654 (1996). [CrossRef]

], finite element modal analysis [16

J. C. Grant, J. C. Beal, and N. J. P. Frenette, “Finite element analysis of the ARROW leaky optical waveguide,” IEEE J. Quantum Electron. 30, 1250–1253 (1994). [CrossRef]

,17

H. P. Uranus, H. J. W. M. Hoekstra, and E. V. Groesen, “Simple high-order Galerkin finite scheme for the investigation of both guided and leaky modes in anisotropic planar waveguides,” Opt. Quantum Electron. 36, 239–257 (2004). [CrossRef]

], and imaginary distance beam propagation method [18

Y. Tsuji and M. Koshiba, “Guided-mode and leaky-mode analysis by imaginary distance beam propagation method based on finite element scheme,” J. Lightwave Technol. 18, 618–623 (2000). [CrossRef]

], have been built up to tackle various leaky waveguides.

Recently, an efficient and accurate numerical method, based on a number of orthogonal basis functions termed as the pseudospectral method [19

J. P. Boyd, “Chebyshev and Fourier Spectral methods,” in Lecture Notes in Engineering, 2nd ed. (Springer Verlag, 2001).

], has been proposed to solve dielectric optical waveguides supporting only the guided modes [20–22

C. C. Huang, C. C. Huang, and J. Y. Yang, “An efficient method for computing optical waveguides with discontinuous refractive index profiles using spectral collocation method with domain decomposition,” J. Lightwave Technol. 21, 2284–2296 (2003). [CrossRef]

]. The purpose of this paper is to extend the ability of the pseudospectral scheme to the ARROW structures with leaky modes. At the present time, to modify the consideration representing outermost subdomains by the Laguerre Gaussian (LG) functions for guided mode problems with the exponential decay feature in our previous work [20–22

], the outermost subdomains occupied by oscillatory leaky waves are represented by the Chebyshev polynomials with Mur’s absorbing boundary condition (ABC) [23

G. Mur, “Absorbing boundary conditions for the finite difference approximation of the time-domain electromagnetic field equations,” IEEE Trans. Electromagn. Compat. 23, 377–382 (1981). [CrossRef]

]. Accordingly, the oscillatory behaviors of leaky waves extending to infinity can be well represented, and the superior performances of pseudospectral scheme for accuracy and efficiency may also be preserved. A number of isotropic and anisotropic ARROW structures are analyzed to verify the performances of the proposed solution method. Other issues of the paper are organized as follows. In Section 2, we formulate the wave equations for an anisotropic medium. The description of the pseudospectral scheme with Mur’s ABC is stated in Section 3. Section 4 simulates several numerical examples to validate the accuracy and efficiency of the present approach by comparing exact data. Finally, Section 5 concludes this work.

2. Formulations of planar waveguides

Considering an anisotropic medium if the principal dielectric axes of the crystal are found to be parallel to the waveguide coordinate system, the dielectric tensor ε̿ is expressed as follows:

\overset{=}{ε} = ε_{0} [\begin{matrix} n_{xx}^{2} & 0 & 0 \\ 0 & n_{yy}^{2} & 0 \\ 0 & 0 & n_{zz}^{2} \end{matrix}],

(1)

where ε ₀ is the dielectric constant in free space, and n_xx , n_yy , and n_zz are the refractive indices along the x, y, and z axes of the waveguide coordinate system, respectively. In this paper, we aim at analyzing leaky modes in planar ARROW of arbitrary index profiles with a diagonal dielectric tensor shown in expression (1). Assuming that ε̿ is only a function of the independent x (i.e., planar waveguides with finite and infinite extension in x and y directions, respectively), for the source-free, time-harmonic of exp(iωt) form, and nonmagnetic media, the modal equations in frequency domain derived from Maxwell’s two curl equations are given as follows:

\frac{\partial}{\partial x} (\frac{\partial E_{y} (x)}{\partial x}) + k_{0}^{2} (n_{yy}^{2} (x) - n_{ff}^{2}) E_{y} (x) = 0,

(2)

for transverse-electric (TE) and that

n_{xx}^{2} (x) \frac{\partial}{\partial x} (\frac{1}{n_{zz}^{2} (x)} \frac{\partial H_{y} (x)}{\partial x}) + k_{0}^{2} (n_{xx}^{2} (x) - n_{eff}^{2}) H_{y} (x) = 0

(3)

for transverse-magnetic (TM) polarizations. Here, k ₀=2π/λ ₀, λ ₀ denotes the wavelength in free space, n_eff denotes the complex effective refractive index, and E_y (x) and H_y (x) represent, respectively, the components of the electric and magnetic fields in the y -direction.

3. Numerical scheme

3.1 Pseudospectral scheme and interfacial patching conditions

The framework of the pseudospectral scheme is to divide computational domain into a number of subdomains determined by discontinuous material interfaces. In each subdomain, the optical field φ(x) is expanded by a set of proper interpolation functions C_k (x) and the unknown grid point values φ_k as follows:

φ (x) = \sum_{k = 0}^{n} C_{k} (x) φ_{k},

(4)

where

C_{k} (x) = \frac{ψ_{n + 1} (x)}{ψ_{n + 1}^{'} (x) (x - x_{k})}, 0 \leq k \leq n .

(5)

Here, ψ _n+1(x)denotes an arbitrary basis function of order n+1 and is determined according to the feature of that subdomain, the prime denotes the first derivative of ψ _n+1(x) with respect to x, and x_k denotes the corresponding collocation point fulfilling the condition of C_l (x_k )=δ_lk , where δ_lk denotes the Kronecker delta. The explicit forms of C_k (x) for distinct basis functions are represented below. If we take the Chebyshev polynomials as basis functions, we have [19

J. P. Boyd, “Chebyshev and Fourier Spectral methods,” in Lecture Notes in Engineering, 2nd ed. (Springer Verlag, 2001).

]

C_{k} (x) = \frac{{(- 1)}^{k + 1} (1 - x^{2}) T_{n}^{'} (x)}{c_{k} n^{2} (x - x_{k})}, x \neq x_{k}

(6)

where T_n (x) denotes the Chebyshev polynomial of order n, c ₀=c_n =2, and c_k =1 (1≤k≤n-1). For the LG functions, namely ψ _n+1(αx)=exp(-αx/2)(αx)L_n (αx) whereL_n (αx)denotes the Laguerre polynomial of order n, the explicit form ofC_k (x)is given as follows [19

J. P. Boyd, “Chebyshev and Fourier Spectral methods,” in Lecture Notes in Engineering, 2nd ed. (Springer Verlag, 2001).

C_{k} (α x) = \frac{e^{- \frac{α x}{2}}}{e^{- \frac{α x_{k}}{2}}} \frac{(α x) L_{n} (α x)}{(α x L_{n})' (α x_{k}) (α x - α x_{k})}, x \neq x_{k} .

(7)

Here, in LG functions, exists an important parameter α called as a scaling factor that can influence the accuracy for a given number of terms of basis functions. The definite procedure for determining α has been proposed in detail in our previous work [20

, 21

C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron. 11, 457–465 (2005). [CrossRef]

]. From Eqs. (2) and (3), we have an eigenvalue problem in the form

A Φ = {(k_{0} n_{eff})}^{2} I Φ,

(8)

where I denotes a unit matrix, Φ denotes the vector of unknown optical field, and A denotes the transverse differential operator. For obtaining the entries of the transverse operator to the e th subdomain, we substitute Eq. (4) with a suitable basis function into Eqs. (2) and (3) and obtain

A_{lk}^{e} = C_{k}^{(2)} (x_{l}) + k_{0}^{2} n_{yy}^{2} (x_{l}) δ_{lk}, (l, k = 0, 1, 2, . . . n)

(9)

for TE and that

A_{lk}^{e} = \frac{n_{xx}^{2} (x_{l})}{n_{zz}^{2} (x_{l})} [C_{k}^{(2)} (x_{l}) - \frac{2}{n_{zz} (x_{l})} C_{k}^{(1)} (x_{l})] + k_{0}^{2} n_{xx}^{2} (x_{l}) δ_{lk}, (l, k = 0, 1, 2, . . . n)

(10)

for TM polarizations. In Eqs. (9) and (10),

C_{k}^{(n)}

(x_l ) means the n th derivative result of C_k (x) with respect to x at collocation point x_l . The global transverse operator A can be obtained by assembling the entries of the transverse operators to all subdomains. However, in matrix A, the rows located in the dielectric interfaces x_r ({r=1, 2,…m}, where m denotes the number of discontinuous nodes) shared by two adjacent subdomains are imposed to be replaced by the interfacial boundary conditions

E_{y} (x_{r}^{+}) = E_{y} (x_{r}^{-}), \frac{\partial E_{y} (x_{r}^{+})}{\partial x} = \frac{\partial E_{y} (x_{r}^{-})}{\partial x}

(11)

for TE and that

H_{y} (x_{l}^{+}) = H_{y} (x_{l}^{-}), \frac{n_{zz}^{2} (x_{r}^{-}) \partial H_{y} (x_{r}^{+})}{\partial x} = \frac{n_{zz}^{2} (x_{r}^{+}) \partial H_{y} (x_{r}^{-})}{\partial x}

(12)

for TM polarizations, where the

x_{r}^{+}

and

x_{r}^{-}

are, respectively, referred to the locations at the infinitesimal right and left of the interface x_r .

3.2 Determining basis functions and boundary conditions

The best choice is that the interior and outermost subdomains are expanded, respectively, by the Chebyshev polynomials and LG functions, if only guided modes exist [20–22

]. However, the outermost subdomains represented by the LG functions are inappropriate while tackling the waveguide structures with leaky modes. In this work, we mainly propose an efficient solution method for handling leaky waveguides. In general, the optical fields in the outermost subdomains exhibit two conditions: one is with evanescent wave, and the other is with traveling wave. Certainly, the outermost subdomains with the evanescent decay characteristic can be efficiently represented by the LG functions as well as those considered in Refs. [20–22

]. As for that with traveling field, since the oscillatory behavior resulting from leaky modes penetrating to substrate, which have a larger refractive index than that in the core layer, this kind of subdomain is expanded by the Chebyshev polynomials at collocation points but has to be fulfilled by the Mur’s ABC [23

] at the boundary point (i. e., outmost collocation point), which allows the fields to radiate freely out of the computational window. The perfectly matched layer (PML) condition may be a more efficient scheme than Mur’s ABC, but the implementation of PML to the present method is very complicated. In this paper, we mainly aim to validate the accuracy and efficiency of the pseudospectral approach to leaky mode problems. As a result, we prefer choosing a simpler boundary condition like Mur’s ABC to accomplish the work.

In this work, we adopt the second-order formula of Mur’s schemes in frequency domain, and it can be given as follows:

{(\frac{\partial}{\partial r} + i k_{r})}^{2} φ = 0,

(13)

where

i = \sqrt{- 1}

, φ denotes the unknown field and is expanded by the Chebyshev polynomials, r denotes the outward direction normal to the boundary, and k_r denotes the wavenumber in r-direction. In the paper, r and k_r can be replaced by x and k_x in the present coordinate system, respectively. Further, Eq. (13) can be expanded as follows:

(\frac{\partial^{2}}{\partial x^{2}} + 2 i k_{x} \frac{\partial}{\partial x} - k_{x}^{2}) {φ ∣}_{x = x_{b p}} = 0,

(14)

where x_bp denotes the boundary point in the outermost subdomain. Note that the chosen computational interval of the outermost subdomain expanded by the Chebyshev polynomials slightly affects the convergence. Observing Eq. (14), we find that the plane wave solution of φ(x)=φ(0)exp(-ik_xx)is fulfilled [17

]. After that, we substitute the plane wave solution into Eqs. (2) and (3) for the outermost subdomain with homogeneous refractive index; we have

k_{x} = k_{0} \sqrt{n_{yy}^{2} - n_{eff}^{2}}

(15)

for TE, and

k_{x} = k_{0} \frac{n_{zz}}{n_{xx}} \sqrt{n_{xx}^{2} - n_{eff}^{2}},

(16)

for TM polarizations. Through Eqs. (2) and (3), the first term in Eq. (14) can be replaced by the terms of -

k_{x}^{2}

φ for both TE and TM polarizations. As a result, Eq. (14) can be further modified as follows:

(2 i k_{x} \frac{\partial}{\partial x} - 2 k_{x}^{2}) {φ ∣}_{x = x_{b p}} = 0

(17)

for TE and TM polarizations. For ensuring only the outgoing plane wave, we must restrict that the real part of k_x is positive. In numerical implementation, in order to keep the unit matrix I in Eq. (8), the term (

k_{0}^{2}

n_{eff}^{2}

)φ is added simultaneously in two sides of Eq. (17). Therefore a matrix eigenvalue problem is obtained by evading a generalized matrix eigenvalue problem. Since the introduction of the Mur’s ABC, the eigenvalue n_eff resides within the global matrix and is required to be solved by iteration scheme. The computational procedure is that we first assign the initial guess of

n_{eff}^{0}

to A, and then a new value of

n_{eff}^{1}

is found by solving a general linear matrix eigenvalue. Subsequently,

n_{eff}^{1}

is used as a new value to calculate the next one

n_{eff}^{2}

. We do not accept the convergent results until the differences of real and imaginary part of n_eff between the adjacent iterations are simultaneously less than the order of 10^-6.

Besides the consideration of choosing basis functions in the present work, the Chebyshev polynomials with Mur’s ABC can also be used to represent the evanescent field having decaying characteristics at the boundary as well as that used in the finite element scheme with Sommerfeld’s ABC [17

]. However, in our numerical experiments concerning computational efficiency (the required number of basis functions), the preferable choice is still the LG functions. Moreover, no supplementary boundary conditions are demanded for evanescent fields [20–22

], and thus the numerical implementation may be significantly simplified.

4. Numerical results and discussion

In order to verify the performances of the present scheme as applied to various ARROW structures, we first analyze in detail the propagation characteristics of an isotropic example, which has been studied by Baba and Kokubun [6

] using the system interference matrix and by Chen et al. [14

] using TMM with APM. Except for the anisotropy, the second example is the same as the first example. Finally, we study the coupling length and field distributions of even and odd modes for an ARROW-based directional coupler device. The accuracy of these examples calculated by the proposed approach is discussed and compared with other published results.

4.1 Isotropic ARROW structure

The schematic diagrams of the geometry and refractive index profile of the ARROW structure are depicted in Fig. 1. The relevant parameters used are those with an operating wavelength of λ=0.6328µm, the refractive indices are n_a=1for air and n_s=3.85 for substrate, and the thicknesses (refractive indices) for the first cladding, second cladding, and core layers are, respectively, d₁=0.142λ (n₁=2.3), d₂=3.15λ (n₂=1.46), and d_c=6.3λ (n_c=1.46) under the antiresonant condition. The computational domain is divided into five subdomains. The calculated results of the complex effective indices for the first i time iterations

n_{eff}^{i}

are shown in Table 1, which also shows that

n_{eff}^{0}

=1.459 denotes the chosen initial guess of effective index. In the example, the number of terms of LG basis functions for the outermost subdomain occupied by air is N_a=10 while the computational interval estimated according to our criterion [21

] for spreading of guided mode is S_a =1µm, the number of Chebyshev polynomials for the first cladding, second cladding, and core layers are N₁=N₂=N_c=20, and the Chebyshev polynomials for substrate layer is N_s=40 while the computational interval of the outermost subdomain chosen is S_sub =1µm. More terms of basis functions in the substrate layer are necessary due to the sharply oscillatory behavior of leaky waves. In addition, considering the effect caused by the computational interval of the outermost subdomain with higher refractive index, the differences of the convergent results between S_sub =1µm and S_sub =2µm for all of the modes are negligible while using N_s=50 for S_sub =2µm. We conclude that the major factor in obtaining accurate results is use of the pseudospectral scheme. Table 1 also shows that with only three iterations, our method achieves the convergent values. In addition, the exact four-digit values for n_eff can be obtained even when only one iteration is executed. In our numerical experiments, the present approach obtains the same convergent rate for using different

n_{eff}^{0}

, as long as the chosen value of

n_{eff}^{0}

is slightly smaller than n_c. The solutions of this work, obtained by the third iteration, and the results calculated by Chen et al. [14

] are shown in Table 2. We can see that the two schemes are in excellent agreement.

Fig. 1. Schematic diagrams of geometry and refractive index profile of the ARROW structure.

Table 1. Complex Effective Indices of the TE and TM Modes of an Isotropic ARROW Structure by the Present Method for the First Three Iterations

(N_s=40, N₁=N₂=N_c=20, N_a=10, $n_{eff}^{0}$ =1.459)
Mode	Re( $n_{eff}^{1}$ )	-Im( $n_{eff}^{1}$ )(×10^-3)	Re( $n_{eff}^{2}$ )	-Im( $n_{eff}^{2}$ )(×10^-3)	Re( $n_{eff}^{3}$ )	-Im( $n_{eff}^{3}$ )(×10^-3)
TE₀	1.802689419	≃0	1.802689419	≃0	1.802689419	≃0
TE₁	1.457941265	0.000054189	1.457941265	0.000054189	1.457941265	0.000054189
TE₂	1.451919137	0.052846592	1.451919174	0.052870681	1.451919174	0.052870681
TE₃	1.451173877	0.192047522	1.451174055	0.192035341	1.451174055	0.192035341
TE₄	1.441371354	0.004372767	1.441371363	0.004374469	1.441371363	0.004374469
TE₅	1.427413547	0.213184218	1.427414119	0.213733418	1.427414119	0.213733398
TE₆	1.424444291	0.766252682	1.424447391	0.766727773	1.424447391	0.766727684
TM₀	1.585174751	0.000000004	1.585174679	0.000000004	1.585174679	0.000000004
TM₁	1.457890856	0.002450742	1.457890856	0.002450742	1.457890856	0.002450742
TM₂	1.451755121	0.552864213	1.451754691	0.553891766	1.451754691	0.553891897
TM₃	1.451302333	1.152697420	1.451304282	1.151032482	1.451304282	1.151033285
TM₄	1.440916234	0.190565224	1.440916633	0.190617140	1.440916633	0.190617138
TM₅	1.426481018	2.169331527	1.426462968	2.181427554	1.426462961	2.181440421
TM₆	1.425511013	5.078997130	1.425568355	5.056058275	1.425568264	5.056121818

The initial guess of effective index is

n_{eff}^{0}

=1.459, and the number of terms of basis function in each subdomain is N_a=10, N₁=N₂=N_c=20, and N_s=40.

Table 2. Calculated Results of the Present Method and TMM Using APM [14

]

This work			TMM using APM [14 C. K. Chen, P. Berini, D. Feng, S. Tanev, and V. P. Tzolov, “Efficient and accurate numerical analysis of multilayer planar waveguides in lossy anisotropic media,” Opt. Express 7, 260–272 (2000). [CrossRef] [PubMed] ]
Mode	Re(n_eff )	-Im(n_eff )(×10^-3)	Re(n_eff )	-Im(n_eff )(×10^-3)
TE₀	1.802689419	≃0	-
TE₁	1.457941265	0.000054189	1.457941265	0.000054189
TE₂	1.451919174	0.052870681	1.451919174	0.052870681
TE₃	1.451174055	0.192035341	1.451174055	0.192035341
TE₄	1.441371363	0.004374469	1.441371363	0.004374469
TE₅	1.427414119	0.213733398	-
TE₆	1.424447391	0.766727684	-
TM₀	1.585174679	0.000000004	-
TM₁	1.457890856	0.002450742	1.457890856	0.002450742
TM₂	1.451754691	0.553891897	1.451754691	0.553891897
TM₃	1.451304282	1.151033285	1.451304282	1.151033285
TM₄	1.440916633	0.190617138	1.440916633	0.19061714?
TM₅	1.426462961	2.181440421	-
TM₆	1.425568264	5.056121818	-

To further demonstrate the accuracy of our method, we calculate the dispersion and radiation loss characteristics of the ARROW structure versus the thickness of each layer while the thicknesses of other layers are at corresponding antiresonant conditions. First, Figs. 2(a) and 2(b) show the dispersion and loss characteristics of the first six ARROW modes of TE polarization, labeled TE₁-TE₆, versus the thickness of the first cladding (d₁/λ), respectively, while other parameters are fixed. Here, referring to Ref. [6

], we call TE₀ as the first cladding mode. In Fig. 2(a), the ranges of d₁/λ with flat portions of effective index for each mode correspond to the antiresonances of Fabry-Perot cavities formed in the two interference claddings. The same portions represent the low-loss portions while referring to Fig. 2(b). It is clear that the low-loss portion of TE1 is broad so as to allow large fabrication tolerance. In contrast to the flat portions, the transitional portions in Fig. 2(a) illustrate the resonances subject to high-loss portions in Fig. 2(b). For comparing the order of losses between different polarized waves, Fig. 2(b) also illustrates the first two ARROW modes of TM polarization. It is clear that the losses of TM modes are larger than TE ones at a corresponding mode number. As a result, the higher-order TE and TM modes can be easily filtered out to obtain single-mode (TE₁) propagation. Besides, in Fig. 2(b), the loss characteristics display periodic variety as a function of (d₁/λ). Likewise, the propagation characteristics versus the thickness of the second cladding (d₂/λ) and core (d_c/λ) layers are shown in Figs. 3 and 4, respectively. In Figs. 3(b) and 4(b), at corresponding first antiresonant conditions, high loss discrimination as that observed in Fig. 2(b) is also shown.

Fig. 2. (a). Dispersion characteristics; (b).Radiation loss characteristics of various modes of an isotropic ARROW versus the thickness of the first cladding.

Fig. 3. (a). Dispersion characteristics; (b). Radiation loss characteristics of various modes of an isotropic ARROW versus the thickness of the second cladding.

Fig. 4. (a). Dispersion characteristics; (b). Radiation loss characteristics of various modes of an isotropic ARROW versus the thickness of the core.

Furthermore, Figs. 5(a)–5(f) show the real part of relative field profiles of various TE (solid curve) and TM (dotted curve) polarized waves simultaneously at the first antiresonant condition (i. e., d₁=0.142λ, d₂=3.15λ, and d_c=6.3λ). Here, “relative” means that the field profiles are normalized by their maximum values [6

]. In this paper, Figs. 5(a)–5(f) are sequentially labeled as TE₀ (TM₀) to TE₆ (TM₆). We can observe that the fields in Fig. 5(a) are confined to the vicinity of the first cladding layer, and the confinement of the TE mode is tighter than that of the TM mode. Figure 5(b) illustrates that the lowest ARROW mode (TE₁) is most localized in the core layer, which is of practical interest, and that it is similar to the fundamental mode of conventional waveguides. The oscillatory characteristics in the substrate layer are also clearly shown, and the larger loss of TM mode is observed, as well as the values listed in Table 2. From Figs. 5(b)–5(f), excluding Figs. 5(b) and 5(e), the fields penetrating the interference claddings are fairly large.

Fig. 5. Relative field profiles of TE and TM modes of isotropic ARROW on the first antiresonant condition on (d₁=0.142λ,d₂=3.15λ, and d_c=6.3λ): (a) TE₀ and TM₀ (the first cladding mode); (b) TE₁ and TM₁; (c)TE₂ and TM₂; (d) TE₃ and TM₃; (e) TE₄ and TM₄; (f) TE₅ and TM₅.

4.2 Anisotropic ARROW structure

Following the same structure analyzed in Subsection 4.1, but with anisotropy in inner layers sandwiched between air and substrate layers, the anisotropic material with a diagonal form of dielectric tensor is expressed in Eq. (1). The refractive indices are n_zzc=1.46 for the core, n_zz1=2.3for the first cladding, and n_zz2=1.46 for the second cladding. Other values follow the relations of n_xxj=n_yyj=1.03n_zzj, where (j=c, 1,2). Here, the initial guess of effective index of

n_{eff}^{0}

=1.459 is chosen, and three iterations are executed to achieve convergent solutions. The results obtained by our method, using 140 unknowns (N_a=10,N₁=N₂=N_c=30,N_s=40), are shown in Table 3 with the six-order accuracy of the finite element method [17

] and TMM using APM [14

]. The number of unknowns in Ref. [17

] is nearly 1300 because of the effective mesh size of 0.005µm. This shows that the computational efficiency of the present scheme is superior to that of the finite element method with six-order accuracy [17

]. In addition, the results of complex effective indices obtained by our work show excellent agreement with exact solutions [14

]. However, the present scheme obtains a full matrix. The computational time is 0.18 sec (each iteration) executed on Pentium IV PC with a CPU clock rate of 3.0 GHz in a double precision. We obtain that the effective indices and losses as a function of thickness of each inner layer are similar to its isotropic counterpart. As for that the principal dielectric axes of the crystal are not parallel to the waveguide coordinate system, the more complicated hybrid modes and the propagation characteristics depend strongly on the included angle between the optics axis and the waveguide coordinate system.

Table 3. Calculated Results of an Anisotropic ARROW Structure by the Present Method, Six-Order Accuracy FEM [17

], and TMM Using APM [14

]

(N_s=40, N₁=N₂=N_c=30, N_a=10, $n_{eff}^{0}$ =1.459)
Mode	This work		Six-order FEM [17 H. P. Uranus, H. J. W. M. Hoekstra, and E. V. Groesen, “Simple high-order Galerkin finite scheme for the investigation of both guided and leaky modes in anisotropic planar waveguides,” Opt. Quantum Electron. 36, 239–257 (2004). [CrossRef] ]		TMM using APM [14 C. K. Chen, P. Berini, D. Feng, S. Tanev, and V. P. Tzolov, “Efficient and accurate numerical analysis of multilayer planar waveguides in lossy anisotropic media,” Opt. Express 7, 260–272 (2000). [CrossRef] [PubMed] ]
Mode	Re(n_eff )	-Im(n_eff )(×10^-3)	Re(n_eff )	-Im(n_eff )(×10^-3)	Re(n_eff )	-Im(n_eff )(×10^-3)
TE₀	1.867833877	≃0	1.867833876	≃0	-
TE₁	1.501798936	0.000050179	1.501798936	0.000050179	1.501798936	0.000050179
TE₂	1.495945499	0.053815143	1.495945499	0.053815130	1.495945499	0.053815143
TE₃	1.495255344	0.184243873	1.495255344	0.184243812	1.495255344	0.184243873
TE₄	1.485698165	0.004051178	1.485698165	0.004051177	1.485698165	0.004051178
TE₅	1.472140044	0.217484450	-		-
TE₆	1.469389780	0.735702355	-		-
TM₀	1.632729920	0.000000004	1.632729919	0.000000004	-
TM₁	1.501625054	0.002544521	1.501625054	0.002544521	1.501625054	0.002544521
TM₂	1.495287895	0.576101022	1.495287895	0.576101022	1.495287895	0.576101022
TM₃	1.494855078	1.189339701	1.494855078	1.189339701	1.494855078	1.189339701
TM₄	1.484121307	0.197863211	1.484121307	0.197863211	1.484121307	0.197863211
TM₅	1.469174801	0.224772521	-		-
TM₆	1.468390303	5.256321451	-		-

4.3 ARROW-based directional coupler

The directional coupler, which consists of two parallel waveguides, is an interesting photonic structure to build various optical communication devices such as a modulator, filter, power divider, and polarizer. The coupling mechanism of a conventional directional coupler is well studied by the coupled mode theory [24

A. Hardy and W. Streifer, “Coupled mode theory of parallel waveguides,” J. Lightwave Technol. 3, 1135–1146 (1985). [CrossRef]

], which considers merely the weakly coupled approximation. In addition, the coupling strength decreases as an exponential function of waveguide separation. Because of the widespread applications of ARROW-based structures, a novel dual ARROW-based directional coupler utilizing antiresonant reflection as the guiding mechanism has been proposed to investigate the propagation characteristics [25

M. Mann, U. Trutschel, C. Wachter, L. Leine, and F. Lederer, “Directional coupler based on an antiresonant reflecting optical waveguide,” Opt. Lett. 16, 805–807 (1991). [CrossRef] [PubMed]

]. In comparison with a conventional directional coupler that utilizes total internal reflection, the ARROW structures have strong coupling, which is due to their having leaky waves rather than the evanescent waves found in conventional two parallel waveguides; furthermore, the variation of the coupling length versus waveguide separation displays periodicity [25

M. Mann, U. Trutschel, C. Wachter, L. Leine, and F. Lederer, “Directional coupler based on an antiresonant reflecting optical waveguide,” Opt. Lett. 16, 805–807 (1991). [CrossRef] [PubMed]

In Fig. 6, the configuration of the dual ARROW structure grown on a Si substrate with a refractive index n_s=3.5 consists of two ARROW structures separated by a separation cladding layer with a refractive index n_sep=1.46 and thickness d_sep=2µm [26