Higher order finite-difference frequency domain analysis of 2-D photonic crystals with curved dielectric interfaces

Yen-Chung Chiang

doi:10.1364/OE.17.003305

1. Introduction

Although the term “photonic crystal” (PC) was first used over one century ago, the idea that we might be able to control the properties of light propagation just as we do in semiconductor devices did not take off until Eli Yablonovitch [1] and Sajeev John [2] published two milestone papers on PCs in 1987. PCs consist of periodic dielectric or metallo-dielectric structures, and they are characterized by their photonic band gaps (PBGs) in which waves are forbidden to propagate in all directions for some frequencies. Based on this characteristic, we may design various devices by properly introducing defects between the periodic structures. Thus the calculation of band diagrams is a fundamental task for analyzing PCs and designing devices.

Among various methods proposed for calculating the band structures of the PCs, the most commonly used ones are the plane-wave expansion (PWE) method [3–5] and the finite-difference time-domain (FDTD) method [6–8]. Aside from these two methods, Yang [9] proposed a finite-difference frequency domain (FDFD) scheme for analyzing the band structures of two dimensional (2-D) PCs, and he also demonstrated the sparsity of the resulting matrix and the efficiency of this method. But such scheme is suitable only for the irregular shapes with interfaces parallel to the axes but not for those with arbitrary geometries. Conventional finite-difference methods generally adopt the staircase-index [9–11], the gradual-index [12], or the index-average [13] approximation to model the slanted or curved interfaces. But these approximations cannot treat interfaces of high index contrast, that is, the strongly guided case very well. Improved finite-difference schemes considering the step-index interface have also been proposed [14–16] to deal with 1-D problems, but these methods fail to treat 2-D structures with the slanted or curved interfaces that are usually found in PC structures. Although some finite-difference schemes have been proposed for the tilt step-index interfaces [17], approximations as the index-average assumption or interpolation in geometry are still required in those schemes. To overcome the difficulty of adopting finite-difference method in structures with curved step-index interfaces, suitable finite difference schemes [18, 19] by a systematic derivation procedure have been proposed for modal analysis and band diagram calculation, respectively.

According to the FDFD idea provided by Yang [9], Yu and Chang [20] recently also proposed a compact FDFD method for analyzing the band structures of PCs based on Yee-mesh formulation and index-average scheme [13]. The results offered by [20] seemed to increase numerical efficiency and accuracy, but it was reported [21] that it suffered from non-uniform numerical convergent speed between different bands. In addition to the non-uniform convergence property, the final formulation in [20] involves global matrix inversion and the resulting matrix will lose its sparsity, which is the most attractive factor of adopting the FDFD method. With the help of the generalized Douglas scheme [16] or the Padé approximation [22], we will extend our previous works in [19] to high-order formulas for analyzing the band structures of PCs. The procedure of deriving high-order finite-difference formulations is given in Section 2. Some numerical results are presented in Section 3 for PCs with rectangular and triangular lattice cases including the comparison of our method with previous works. The conclusion is drawn in Section 4.

Fig. 1. The cross-sectional view of a 2-D PC and its unit cell with a being the lattice distance and r being the radius of the circular rods for (a) the square lattice and (b) the triangular lattice.

Download Full Size | PDF

Fig. 2. Cross section of general finite-difference mesh with a curved interface.

Download Full Size | PDF

2. The high-order FDFD scheme

The cross-section of a 2-D PC is shown in Figs. 1(a) and (b) for the square lattice and the triangular lattice cases, respectively. The unit cells for both cases are also illustrated as the regions enclosed by blue dashed lines in Figs. 1(a) and (b), where a and r denote the lattice distance and the radius of element rods, respectively. The element rods can be either dielectric rods or air columns. As indicated in [9], the propagation in the PCs can be categorized into in-plane and out-of-plane cases. Since the 2-D PC is uniform along the z direction and periodic in the transverse plane, we will only consider the in-plane propagation for calculating the band structures and let the propagation constant in the z direction, that is, β_z be zero.

2.1. High-order relations between grid points

Consider the cross section shown in Fig. 2, in which a curved interface lies between two dielectric regions with permittivities of ε_L and ε_R, respectively. The quantity ϕ denotes E_z for TM modes and H_z for TE modes at grid points. A local cylindrical coordinate system with origin O′ can be defined in terms of the normal vector r̂, the tangential vector θ̂, and the effective radius R. To derive a correct high-order finite-difference formula, the general relation between the field at a sampled point, such as ϕ|_(m,n) in Fig. 2, and the fields at nearby points, must first be found. The basic procedure for deriving this relation is similar to those provided in the previous work [19]. Using the field relation between ϕ|_(m+1,n+1) and ϕ|_(m,n) as an example, the steps that link the fields at two grid points can be illustrated as

ϕ ∣_{(m + 1, n + 1)} \overset{TSE}{\leftrightarrow} ϕ ∣_{R; x - y} \overset{LCT}{\leftrightarrow} ϕ ∣_{R; r - θ} \overset{BCM}{\leftrightarrow} ϕ ∣_{L; r - θ} \overset{LCT}{\leftrightarrow} ϕ ∣_{L; x - y} \overset{TSE}{\leftrightarrow} ϕ ∣_{(m, n)},

where TSE, LCT, and BCM denote the 2-D Taylor series expansion, local coordinate transformation, and boundary conditions matching, respectively. Since the manipulation of the Taylor series expansion and local coordinate transformation are relatively easy and the lower order equations can be found in [19], the related high-order equations for TSE and LCT are omitted to save space.

As for the boundary conditions matching step, the quantity ϕ|_R;r-θ, which is the field just to the right of the interface, and its derivatives must be linked with those on the other side of the interface. The required equations with derivative terms up to second order can also be found explicitly in [19]. To derive other high-order terms, we notice that the interface is laid in the tangential direction in the local coordinate system, as shown in Fig. 2. If a function is continuous along the tangential variable θ, then its derivatives with respect to θ must also be continuous in this tangential direction. Therefore, if a low-order boundary matching equation of the form

\frac{\partial^{p} ϕ}{\partial r^{q} \partial θ^{p - q}} ∣_{R} = \tilde{ϕ} (r, θ) ∣_{L}

is available, then the following relation, obtained by differentiating Eq. (2) with respect to θ will also apply.

\frac{\partial^{p + 1} α}{\partial r^{q} \partial θ^{p - q + 1}} ∣_{R} = \frac{\partial \tilde{ϕ} (r, θ)}{\partial θ} ∣_{L}

From the second-order boundary matching equations provided in [19] and the relation between Eqs. (2) and (3), we can obtain the following equations

\frac{\partial^{3} E_{z}}{\partial r^{2} \partial θ} ∣_{R} = \frac{\partial^{3} E_{z}}{\partial r^{2} \partial θ} ∣_{L} + k_{0}^{2} (ε_{L} - ε_{R}) \frac{\partial E_{z}}{\partial θ} ∣_{L}

\frac{\partial^{3} E_{z}}{\partial r \partial θ^{2}} ∣_{R} = \frac{\partial^{3} E_{z}}{\partial r \partial θ^{2}} ∣_{L}

\frac{\partial^{3} E_{z}}{\partial θ^{3}} ∣_{R} = \frac{\partial^{3} E_{z}}{\partial θ^{3}} ∣_{L}

for the TM modes, and

\frac{\partial^{3} H_{z}}{\partial r^{2} \partial θ} ∣_{R} = \frac{\partial^{3} H_{z}}{\partial r^{2} \partial θ} ∣_{L} + \frac{1}{r} (1 - \frac{ε_{R}}{ε_{L}}) \frac{\partial H_{z}}{\partial r} ∣_{L} + k_{0}^{2} (ε_{L} - ε_{R}) \frac{\partial H_{z}}{\partial θ} ∣_{L}

\frac{\partial^{3} H_{z}}{\partial r \partial θ^{2}} ∣_{R} = \frac{ε_{R}}{ε_{L}} \frac{\partial^{3} H_{z}}{\partial r \partial θ^{2}} ∣_{L}

\frac{\partial^{3} H_{z}}{\partial θ^{3}} ∣_{R} = \frac{\partial^{3} H_{z}}{\partial θ^{3}} ∣_{L}

for the TE modes, in which k ₀ is the free-space wavenumber. Besides of the above relations, high-order derivative equations with respect to r only, that is, ∂^pϕ/∂r^p must still be obtained. They can be found by combining the Helmholtz equation with Maxwell’s equations or the Helmholtz equation itself. For example, combining the Helmholtz equation on E_z with the curl equation ∇ × E = - jωμ H yields

\frac{\partial^{3} E_{z}}{\partial r^{3}} ∣_{R} = \frac{\partial^{3} E_{z}}{\partial r^{3}} ∣_{L} + k_{0}^{2} (ε_{L} - ε_{R}) \frac{\partial E_{z}}{\partial r} ∣_{L} - \frac{1}{r} k_{0}^{2} (ε_{L} - ε_{R}) E_{z} ∣_{L} .

Similarly, combining the Helmholtz equation on H_z with another curl equation ∇ × H = jωε E yields

\frac{\partial^{3} H_{z}}{\partial r^{3}} ∣_{R} = \frac{ε_{R}}{ε_{L}} {\frac{\partial^{3} H_{z}}{\partial r^{3}} ∣_{L} + \frac{1}{r} (1 - \frac{ε_{L}}{ε_{R}}) \frac{\partial^{2} H_{z}}{\partial r^{2}} ∣_{L} + [k_{0}^{2} (ε_{L} - ε_{R}) + \frac{1}{r^{2}} (1 - \frac{ε_{L}}{ε_{R}})] \frac{\partial H_{z}}{\partial r} ∣_{L}}

Other high-order relations can also be found by successively applying the above steps. Therefore the relations required for each part in Eq. (1) can be obtained and extended to high-order terms.

2.2. High-order formulation with the GD scheme

Following the steps that are summarized by Eq. (1), the field quantity at each grid point is expressed as a linear combination of the field quantity at point (m,n), ϕ|_(m,n), and its derivatives. Similar linear equations at N grid points, including the central sampled point (m, n) itself, can be expressed in matrix form as,

Φ = M \cdot D_{(m, n)} + H . O . T .,

where Φ is the vector that contains the fields at N grid points; M is an N-by-N matrix that contains coefficients obtained using the steps described in the preceding subsection, and D _(m,n) is the vector that contains the field quantity at point (m, n) and its derivatives with respect to x or y, respectively, up to N required terms. Unlike in Eq. (28) in [19], the high-order terms (H.O.T.) in Eq. (12) are left for further derivation. A set of finite difference formulas are obtained by taking the inverse operation of Eq. (12):

D_{(m, n)} ≅ M^{- 1} \cdot Φ .

If the dimension of the matrix M is 9-by-9, then the nine-point second-order formulation as provided in [19] is obtained. A simple means of obtaining a high-order formula involes the use of more sampled points and more derivative terms in Eq. (13), as adopted in [19] for the 21-point formula. However, the inversion of a matrix for larger dimension in Eq. (13) reduces the efficiency of this method and programming for the case of a PC with the triangular lattice is tedious because of the hexagonal periodic conditions. Therefore, a high-order formulation using only nine sampled points increases the efficiency and reduces the memory requirement.

Notably, the high-order residual terms of the finite-difference scheme based on Eq. (13) can be obtained by multiplying the resultant formulation with Eq. (12) and extracting the terms related to H.O.T. part. Then, the GD scheme introduced in [16, 22] is utilized for treating the residual high-order terms. However, the final form of the fourth-order formula suggested in these previous works is,

\frac{\partial^{2}}{\partial x^{2}} = \frac{D_{xx}^{″}}{1 + c_{1} D_{x}^{'} + c_{2} D_{xx}^{″}},

where D_x′ and D_xx″ denote the first and the second-order finite-difference formulas for the differentiation with respect to x and are given by Eq. (13). This form implies that differentiations with respect to x and y can be decoupled. However, this implication does not apply when a slanted or curved step-index interface lies between grid points. To derive a correct formula, both second-order terms differentiated with respect to x and y must be taken into account simultaneously. Use the field quantity E_z as an example and assume that the formula obtained from Eq. (13) can be expressed as following with the residual third-order terms

(D_{xx}^{″} + D_{yy}^{″}) E_{z} = \frac{\partial^{2} E_{z}}{\partial x^{2}} + \frac{\partial^{2} E_{z}}{\partial y^{2}} + g_{1} \frac{\partial^{3} E_{z}}{\partial x^{3}} + g_{2} \frac{\partial^{3} E_{z}}{\partial y^{3}},

in which the coefficients g ₁ and g ₂ may come from both D_xx″ and D_yy″. Although ∂³ E_z/∂x ²∂y and ∂³ E_z/∂x∂y ² seem to be third-order terms, it should be included in the nine-point second-order formulas using Eq. (13). Thus D_xxy‴ and D_xyy‴ are elements of D _(m,n) and their second-order formulas can be obtained from Eq. (13). We then add proper third-order terms at both sides of Eq. (15) as

(D_{xx}^{″} + D_{yy}^{″}) E_{z} + g_{1} \frac{\partial^{3} E_{z}}{\partial x \partial y^{2}} + g_{2} \frac{\partial^{3} E_{z}}{\partial x^{2} \partial y} = (1 + g_{1} \frac{\partial}{\partial x} + g_{2} \frac{\partial}{\partial y}) (\frac{\partial^{2} E_{z}}{\partial x^{2}} + \frac{\partial^{2} E_{z}}{\partial y^{2}}) .

Rearrange Eq. (16) and replace the partial differential operator by the finite difference formulas obtained from Eq. (13), a final expression will be reached as

(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) E_{z} = (\frac{D_{xx}^{″} + D_{yy}^{″} + g_{1} D_{xyy}^{‴} + g_{2} D_{xxy}^{‴}}{1 + g_{1} D_{x}^{'} + g_{2} D_{y}^{'}}) E_{z} .

Similar derivation process with additional fourth-order residual terms will lead to a fourth-order formula using only nine sampled points.

2.3. Implementation of the high-order formulation

Because of the periodicity in the geometry, we only need to calculate one unit cell area along with proper periodic boundary conditions (PBCs). The unit cells for the rectangular lattice and the triangular lattice are re-drawn with the periodic boundary conditions in Fig. 3(a) and (b). The PBCs for the PC with square lattice can be expressed as

PBC-X : ϕ (x + a, y) = \exp (- j k_{x} a) ϕ (x, y)

PBC-Y : ϕ (x, y + a) = \exp (- j k_{y} a) ϕ (x, y)

where k_x and k_y are the wavenumbers in the x and y directions, respectively.

For the PC with triangular lattice shown in Fig. 3(b), we have the PBCs

PBC 1 : ϕ (x + \frac{\sqrt{3} a}{2}, y - \frac{a}{2}) = \exp [- j (k_{x} \frac{\sqrt{3} a}{2} - k_{y} \frac{a}{2})] ϕ (x, y)

PBC 2 : ϕ (x + \frac{\sqrt{3} a}{2}, y + \frac{a}{2}) = \exp [- j (k_{x} \frac{\sqrt{3} a}{2} + k_{y} \frac{a}{2})] ϕ (x, y)

PBC 3 : ϕ (x, y + a) = \exp (- j k_{y} a) ϕ (x, y) .

Fig. 3. The unit cell of the PC with corresponding PBCs for (a) the square lattice case, (b) the triangular lattice case, and (c) the modified unit cell of (b).

Download Full Size | PDF

The unit cell for the triangular lattice case is hexagonal in geometry, and it is not convenient to be fitted for the general finite difference scheme. In our calculation, we adopted the modified unit cell suggested by [20] as shown in Fig. 3(c), in which the upperright triangle (the magenta one) is shifted to the lowerleft and the lowerright triangle (the cyan one) is shifted to the upperleft. The corresponding boundaries for the modified unit cell still obey the same PBCs as in Eqs. (19a)–(19c).

Applying the high-order finite-difference formulas in the Helmholtz’s equation for all mesh points in the unit cell with the corresponding PBCs, a final equation in the matrix form can be obtained as

[A - {(\frac{ω}{c})}^{2} B] \cdot \tilde{Φ} = 0

where ω is the angular frequency, c is the speed of light, $\tilde{Φ}$ is the global vector whose elements are the fields on all the grid points, A is the characteristic matrix with entries related to the coefficients in the numerator part of Eq. (17), and B is no longer an identical matrix as in Eq. (32) of [19] but with entries related to the coefficients in the denominator part of Eq. (17). The eigenvalues can be obtained by searching the values of ω that makes the matrix determinant zero.

3. Numerical results

The first investigated case is the square-arranged 2-D PC as shown in Fig. 1(a), which consists of parallel alumina rods in the air. The relative permittivity of the alumina is ε _A1 = 8.9 and that of the air is ε _air = 1. The radius of the alumina rod is r = 0.2a, where a is the lattice constant as shown in Fig. 1(a). Each point along the boundary of the first Brillouin zone shown as the middle inset in Fig. 4(a) determines the values of k_x and k_y that we use in Eqs. (18a), (18b) and thus in Eq. (20). The calculated band diagrams of the TE and TM modes are plotted in Figs. 4(a) and (b), respectively. The band structures marked with red circles are the results obtained by using our nine-point high-order finite-difference scheme with 30×30 grid points, and the solid blue lines are the results obtained using the MIT Photonic-Bands package [23] based on the PWE method with 128 × 128 resolution. It can be seen that they match each other quite well for both TE and TM polarizations. The computation is performed on a Pentium IV 1.7 GHz personal computer for our case, and the required CPU time is just a few minutes which is close to that used by adopting nine-point second-order sheme in [19].

To examine the efficiency and the convergence behavior of our high order scheme, we fix the wave vector k at the M point in the first Brillouin zone as the middle inset in Fig. 4(a), i.e., k = (π/a, π/a, 0), and search for the eigen frequencies using different numbers of grid points. We calculated the results of the four most fundamental bands using the nine-point second-order scheme of [19] and the nine-point high-order scheme derived previously as shown in Figs. 5(a)–(d), and the convergence behaviors of results extracted from [20] are also plotted in each mode for comparison. The magenta dash-dot lines in Figs. 5(a)–(d) are the results of average-index method [13] and the red dashed lines are the results using the compact FDFD method [20]. The black dashed lines marked with circles are the results using the nine-point second-order scheme [19], while the solid blue lines are the results using our nine-point high-order scheme proveded here. The relative error in Figs. 5(a)–(d) refers to the percentage of the difference between the calculated results and a reference result obtained with the number of grid points being 2500. This definition is in consistent with those used in [19, 20]. Compared with the results using other methods, our method gives better convergence property. It shows in Figs. 5(a)–(d) that the convergency order for the index-average method in [13] and the nine-point scheme in [19] are of second order. It can also be seen that the non-uniform numerical convergent speed between different bands using the method in [20]. It shows that the results for all four bands using our nine-point high-order scheme can achieve a relative error less than about 0.25%, as the number the of total grid points is larger than 500. To achieve the same accuracy, over 1100 grid points need to be used for the nine-point second-order scheme in [19]. In such case, our high-order formulation reduces the CPU time and memory requirement to about 46% and 42%, respectively, of those used by adopting the nine-point second-order scheme.

Fig. 4. Calculated band diagrams for the 2-D PC formed by square-arranged alumina rods with r/a = 0.2 and ε _A1 = 8.9 in the air. The results are compared with those from PWE method [23]. (a) TE mode and (b) TM mode.

Download Full Size | PDF

Fig. 5. The convergence properties of our method for the 2-D square-lattice PC compared with other methods. (a) TE first band; (b) TE second band; (c) TM first band; (d) TM second band.

Download Full Size | PDF

Fig. 6. Calculated band diagrams for the 2-D PC formed by triangle-arranged dielectric rods with r/a = 0.2 and ε = 11.4 in the air. The results are compared with those from PWE method [23]. (a) TE mode and (b) TM mode.

Download Full Size | PDF

The second investigated case is the triangle-arranged 2-D PC as shown in Fig. 1(b), which consists of parallel dielectric rods with relative permittivity ε = 11.4 and radius r = 0.2a in the air, where a is the lattice constant as shown in Fig. 1(b). Each point along the boundary of the first Brillouin zone shown as the middle inset in Fig. 6(a) determines the values of k_x and k_y that we use in Eqs. (19a)–(19c) and thus in Eq. (20). The calculated band diagrams of the TE and TM modes are plotted in Figs. 6(a) and (b). The band structures marked with red circles are the results obtained using our nine-point high-order finite-difference scheme with 30 × 30 grid points and adopting the modified unit cell as shown in Fig. 3(c). The solid blue lines are the results obtained using the MIT Photonic-Bands package [23]. They again match each other quite well for both TE and TM polarizations. To examine the convergence characteristics, we now consider the wave vector k at the K point in the first Brillouin zone as the middle inset in Fig. 6(a), i.e., k = [2π/(√3a),2π/(3a),0], and use similar processes to calculate the band structure as in the square lattice case. The convergence behaviors of the four most fundamental bands are illustrated in Figs. 7(a)–(d). The lines notations used here are the same as those in Figs. 5(a)–(d). It shows again that our method provides better numerical convergence compared with other methods.

Fig. 7. The convergence properties of our method for the 2-D square-lattice PC compared with other methods. (a) TE first band; (b) TE second band; (c) TM first band; (d) TM second band.

Download Full Size | PDF

Fig. 8. Resulting matrices utilizing this high-order scheme for the calculation of band diagrams of PCs with (a) the rectangular lattice and (b) the triangular lattice.

Download Full Size | PDF

Figs. 8(a) and (b) show the distribution of non-zero elements within the resulting matrices for the square and the triangular lattice cases, respectively. The matrices are obtained by using 30×30 grid points in the unit cell. Since the matrix inversion is only performed on the local matrix in the formula derivation process but not on the global matrix, the resulting matrix remains its sparsity. It can be seen that the matrix is almost band-diagonal with corner elements due to the PBC for the square lattice case. As for the triangle-arranged PC, two side bands appear between diagonal and corner because of the modification in the unit cell as shown in Fig. 3(c). However, the positions of the additional non-zero elements can be determined in advance, thus proper memory assignment and careful programming can still maintain the sparsity of the matrix and reduce the memory requirement and the computational time.

4. Conclusion

We have demonstrated a high-order finite-difference scheme for solving the band diagrams of 2-D PCs with curved interfaces. With the help of the GD scheme, we extended the nine-point second-order shceme to fourth-order accuracy. Our results can match very well with those obtained by PWE method, and it is shown that our scheme can improve convergence behavior significantly. Since we use the local coordinate transformation and boundary conditions matching in the formulation, our higher order scheme is still valid even when slanted or curved step-index interfaces exist in the structure.

Acknowledgments

This work was supported by the National Science Council of the Republic of China under Grant NSC 95-2221-E-005-127 and NSC 97-2221-E-005-091-MY2. It was also supported in part by the Ministry of Education, Taiwan, R.O.C. under the ATU plan. Ted Knoy is appreciated for his editorial assistance.

References and links

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

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

3. K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic dielectric structures,” Phys. Rev. Lett. 65, 3152–3155 (1990). [CrossRef] [PubMed]

4. M. Plihal and A. A. Maradudin, “Photonic band structure of two-dimensional systems: The triangular lattice,” Phys. Rev. B 44, 8565–8571 (1991). [CrossRef]

5. R. D. Meade, K. D. Brommer, A. M. Rappe, and J. D. Joannopoulos, “Existence of a photonic band gap in two dimensions,” Appl. Phys. Lett. 61, 495–497 (1992). [CrossRef]

6. C. T. Chan, Q. L. Yu, and K. M. Ho, “Order-N spectral method for electromagnetic waves,” Phys. Rev. B 51, 16635–16642 (1995). [CrossRef]

7. M. Qiu and S. He, “A nonorthogonal finite-difference time-domain method for computing the band structure of a two-dimensional photonic crystal with dielectric and metallic inclusions,” Appl. Phys. 87, 8268–8275 (2000). [CrossRef]

8. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd edition, (Artech House, MA USA, 2005).

9. H. Y. D. Yang, “Finite difference analysis of 2-D photonic crystals,” IEEE Trans. Microwave Theory Tech. 44, 2688–2695 (1996). [CrossRef]

10. K. Bierwirth, N. Schulz, and F. Arndt, “Finite-difference analysis of rectangular dielectric waveguide structures,” IEEE Trans. Microwave Theory Tech. 34, 1104–1114 (1986). [CrossRef]

11. G. R. Hadley and R. E. Smith, “Full-vector waveguide modeling using an iterative finite-difference method with transparent boundary conditions,” J. Lightw. Technol. 13, 465–469 (1995). [CrossRef]

12. W. P. Huang, C. L. Xu, S. T. Chu, and S. K. Chaudhuri, “The finite-difference vector beam propagation method: analysis and assessment,” J. Lightw. Technol. 10, 295–305 (1992). [CrossRef]

13. Z. Zhu and T. G. Brown, “Full-vectorial finite-difference analysis of microstructured optical fibers,” Opt. Express 10, 853–864 (2002), http://www.opticsexpress.org/abstract.cfm?uri=OE-10-17-853. [PubMed]

14. M. S. Stern, “Semivectorial polarized finite difference method for optical waveguides with arbitrary index profiles,” Proc. Inst. Elect. Eng. J. 135, 56–63 (1988).

15. C. Vassallo, “Improvement of finite difference methods for step-index optical waveguides,” Proc. Inst. Elect. Eng. J. 139, 137–142 (1992).

16. Y.-P. Chiou, Y.-C. Chiang, and H.-C. Chang, “Improved three-point formulas considering the interface conditions in the finite-difference analysis of step-index optical devices,” J. Lightw. Technol. 18, 243–251 (2000). [CrossRef]

17. J. Xia and J. Yu, “New finite-difference scheme for simulations of step-index waveguides with tilt interfaces,” IEEE Photon. Technol. Lett. 15, 1237–1239 (2003). [CrossRef]

18. Y.-C. Chiang, Y.-P. Chiou, and H.-C. Chang, “Improved full-vectorial finite-difference mode solver for optical waveguides with step-index profiles,” J. Lightw. Technol. 20, 1609–1618 (2002). [CrossRef]

19. Y.-C. Chiang, Y.-P. Chiou, and H.-C. Chang, “Finite-difference frequency-domain analysis of 2-D photonic crystals with curved dielectric interfaces,” J. Lightw. Technol. 26, 971–976 (2008). [CrossRef]

20. C.-P. Yu and H.-C. Chang, “Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals,” Opt. Express 12, 1397–1408 (2004), http://www.opticsexpress.org/abstract.cfm?uri=OE-12-7-1397. [CrossRef] [PubMed]

21. P.-J. Chiang, C.-P. Yu, and H.-C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E 75, 026703 (2007). [CrossRef]

22. I. Harari and E. Turkel, “Accurate finite difference methods for time-harmonic wave propagation,” J. Comput. Phys. 119, 252–270 (1995). [CrossRef]

23. S. G. Johnson and J. D. Joannopoulos, “The MIT Photonic-Bands Package home page [on line],” http://ab-initio.mit.edu/mpb/.

Higher order finite-difference frequency domain analysis of 2-D photonic crystals with curved dielectric interfaces

Abstract

1. Introduction

2. The high-order FDFD scheme

2.1. High-order relations between grid points

2.2. High-order formulation with the GD scheme

2.3. Implementation of the high-order formulation

3. Numerical results

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (8)

Equations (24)

Optics Express