1. Introduction

In recent years, photonic crystals (PCs) have inspired a lot of interest and many research efforts have been devoted to their possible applications. They are characterized by their photonic band gaps (PBGs) which result from structure periodicity [1–8]. Two-dimensional (2-D) PCs, such as those composed of either dielectric rods or air columns, are widely employed in many applications such as new means of light waveguiding which offer a very high degree of control over the propagation of light. Even through sharp bends, a wide transmission band can be achieved by carefully designing the geometry of the waveguide. A variety of numerical methods have been utilized to calculate the band structures of 2-D PCs. Among these the most used are the planewave expansion (PWE) method [3–6] and the finite-difference time-domain (FDTD) method [7,8].

The finite difference method (FDM) has been widely adopted to deal with differential equations in numerical analysis. Due to the simplicity of its implementation and the sparsity of its resultant matrix, mode solvers based on the FDM have become attractive numerical methods for analyzing the propagation characteristics of the optical or dielectric waveguides, especially for those complex dielectric waveguides without analytical solutions. The FDM was first employed to solve the scalar waveguide modes under the weakly guiding approximation [9]. For strongly guiding structures, semi-vectorial equations for optical waveguides with arbitrary index profiles were derived [10,11]. To obtain more accurate mode fields, a full-vectorial finite difference scheme was then proposed [12–14]. Each mesh point is surrounded by four sub-regions whose permittivities are allowed to be different. The sparse matrix is formed by discretizing the Helmholtz equations and matching the boundary conditions for all the subregions. Another way is to discretize the wave equations derived from Maxwell’s equations with the finite difference scheme [15–17], which is usually used in the beam propagation method (BPM) [16,17]. Recently, full-vectorial finite-difference waveguide mode solvers based on the Yee’s cell [18] were proposed. Ando et al. [19] derived a Yee-mesh-based imaginary-distance propagation method with modified FD formulae which takes the boundary conditions into account. Fast convergence was achieved for the analysis of step-index fibers and rib waveguides. Zhu and Brown [20] proposed a finite-difference frequency-domain (FDFD) method which directly discretizes Maxwell’s equations using the central difference scheme based on a Yee’s cell and the index average technique was adopted to improve numerical convergence and accuracy.

Besides the PWE method and the FDTD method, finite-difference eigenvalue algorithms for analyzing the band structures of 2-D PCs have been proposed. Yang [21] proposed an algorithm by discretizing the Helmholtz equation in terms of the magnetic field H _z or the electric field E _z for the transverse-electric (TE) and transverse-magnetic (TM) modes, using the central difference scheme in a similar way as that of [12] for waveguide analysis. Only square-lattice PCs were discussed and no curved dielectric interfaces were considered in the unit cell of the PC. Shen et al. [22] derived another algorithm from wave equations in terms of the parallel and transverse components of the magnetic field for the TE and TM modes, respectively, and employed an effective-medium technique for smoothing the field distribution across the dielectric interface. Again, only square PCs were considered. Fast convergent results for 2-D PCs involving dielectric rods were demonstrated.

In this paper, an FDFD eigenvalue algorithm based on the Yee’s cell is derived in the similar way as in [20], and is utilized to analyze the band structures of 2-D PCs with either square or triangular lattice. Instead of discretizing the Helmholtz equations or the wave equations, we treat Maxwell’s equations directly. To increase the accuracy of our FDFD algorithm, a formula based on a fourth-order accurate compact scheme [23] is adopted, and the index average scheme is employed to deal with the curved interface in the computing domain. The derivation of the formulation of this FDFD eigenvalue algorithm is described in Section 2 in detail, and some numerical results are presented and compared with other numerical methods in Section 3. We also discuss the convergence behaviors of our method and other FD methods in Section 4 which is followed by the conclusion.

 

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 circles for (a) the square lattice and (b) the triangular lattice.

Download Full Size | PDF

2. The FDFD algorithm with compact 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, with the respective unit cell illustrated by dashed lines, where a is the lattice distance and r is the radius of the circles which can be either dielectric rods or air columns. Since the 2-D PC is uniform along the z direction and periodic in the transverse plane, in the following discussion we only consider the in-plane propagation that has zero propagation constant in the z direction, β_z=0, so that the wave modes in the PC are either TE or TM to z mode. Assuming all the fields have the exp(jωt) dependence with ω being the angular frequency, Maxwell’s curl equations can be expressed as

where µ and ε are the relative permeability and permittivity of the medium, respectively, and µ ₀ and ε ₀ are the permeability and permittivity of free space.

2.1. The TE mode

For the TE mode, we have only H_z, E_x, and E_y components and the above curl equations can be written as

 

Fig. 2. Yee’s mesh for the (a) TE and (b) TM modes.

Download Full Size | PDF

The transverse plane is then discretized using Yee’s mesh for the TE mode as shown in Fig. 2(a), and we can obtain

where µ_z, ε_x, and ε_y are the relative permeability or permittivity at the points of H_z, E_x, and E_y.

Instead of using the central difference scheme, we introduce a fourth-order compact scheme [23] to replace the differential operator in the above equations to improve the accuracy. For example, $\frac{\partial H_z}{\partial x}$ at points $(i+1,j+\frac{1}{2})$, $(i,j+\frac{1}{2})$, and $(i-1,j+\frac{1}{2})$ can be expanded using the Taylor series expansion with respect to point $(i+\frac{1}{2},j+\frac{1}{2})$ as

where Δx is the grid size in the x direction and H.O.T. means higher order terms. Also, H_z at point $(i-\frac{1}{2},j+\frac{1}{2})$ can be expressed as

If we neglect the higher order terms (H.O.T.) in Eqs. (5)–(8), we can obtain the following relationship between the ∂H_z/∂ x and H_z values

Thus, the term $(\partial H_z⁄\partial x){\mid }_{i,j+\frac{1}{2}}$ in Eq. (4c) can be approximated in terms of H_z at the grid points in the computing domain. Similarly, derivations for ∂H_z/∂ _y, ∂ E_x/∂ _y, and ∂ E_y/∂ _x lead to

With Eqs. (9)–(12), within the computing domain we can obtain the matrix equations

where H_z, E_x, and E_y are the vectors composed of the field components H_z, E_z, and E_y, respectively, at the grid points, and [∂H_z/∂ x], [∂H_z/∂ y], [∂ E_x/∂ y], and [∂ E_y/∂ x] represent the vectors containing the values ∂H_z/∂ x, ∂H_z/∂ y, ∂ E_x∂ y, and ∂ E_y/∂ x, respectively, at the corresponding grid points. U, V, U_x, U_y, V_x, and V_y are square matrices determined by the above compact FDFD formulae and the boundary conditions.

Maxwell’s equations can then be rewritten in a matrix form as

where µ _z, ε _x, and ε _y are diagonal matrices representing µ_z, ε_x, and ε_y at the corresponding grid points. After some mathematical work, we can obtain an eigenvalue equation in terms of H_z as

where k ₀=ω(µ ₀ ε ₀)^1/2 is the wavenumber in free space.

2.2. The TM mode

For the TM mode, we have only E_z, H_x, and H_y components and Maxwell’s curl equations become

which, following the similar derivation using Yee’s cell and the compact scheme, can be converted into the following matrix equation

where the vectors E_z, H_x, and H_y and the diagonal matrices ε _z, µ _x, and µ _y have the similar definitions as in Eq. (13). The eigenvalue equation in terms of E_z for the TM mode can similarly be obtained as

Note that Eqs. (15) and (18) apply to both dielectric and magnetic materials. Therefore, the proposed algorithm is also suitable for calculating the band structures for 2-D PCs composed of magnetic materials.

2.3. Periodic boundary conditions for the 2-D PC

Due to its periodic geometry, in the analysis of the PCs, we only need to consider the unit cell along with the periodic boundary conditions (PBCs). For the PC with square lattice shown in Fig. 1(a), the PBCs can be expressed as

where Ψ is either the electric or magnetic field in the unit cell, and k_x and k_y are the wavenumbers in the x and y directions, respectively

 

Fig. 3. (a) The unit cell of the PC with triangular lattice and its corresponding PBCs. (b) The modified unit cell.

Download Full Size | PDF

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

Please note that Eqs. (20a) and (20b) are sufficient for describing the BCs in an infinitely periodic lattice and Eq. (20c) can be derived from Eqs. (20a) and (20b). PBC1, PBC2, and PBC3 are indicated at the corresponding sides of the unit cell as shown in Fig. 3(a). Since the unit cell for the triangular lattice is hexagonal, it is not convenient to be fitted by Yee’s rectangular mesh. To overcome this problem, an FDTD scheme in a nonorthogonal coordinate system was proposed in [8]. By using the nonorthogonal mesh, the non-rectangular unit cell can be matched quite well but a coordinate transformation is needed and more complicated formulae need to be derived. In our method, we simply move the upperright triangle (the red one) to the lowerleft and the lowerright triangle (the blue one) to the upperleft, as shown in Fig. 3(b), generating a modified unit cell. The modified unit cell still obeys the same PBCs at its boundaries as indicated in Fig. 3(b), e.g., the upper portion of the left side and the lower portion of the right side, both labled as PBC1, are related by the PBC of Eq. (20a). Since it occupies a rectangular area, it is easier to be fitted by the rectangular mesh.

2.4. The average index scheme

As for the curved dielectric interface in the computing domain, the conventional treatment in the FDM is to use either the gradded-index [15–17] or the staircase index [12–14] approximation, and much finer mesh grids are needed to approach the original structure. Improved formulations for oblique or curved interfaces were derived in [24,25] using the Taylor series expansion and matching the interface conditions. Here, for simplicity, an averaging scheme which has been shown to be useful and efficient in increasing the analysis accuracy [20] is adopted to determine the values of µ and ε at each grid point. For example, in Fig. 2(b), ε_z(i, j) is defined as ε_z(i, j)=f·ε ₁ +(1-f)ε ₂, where ε ₁ and ε ₂ are the relative permittivities of media 1 and 2, respectively, in the mesh and f is the filling fraction of medium 1 in the mesh.

3. Numerical results

We first consider the square-arranged 2-D PC, shown as the inset in the middle of Fig. 4(a), formed by parallel alumina rods with relative permittivity ε=8.9 and radius r=0.2a in the air, where a is the lattice distance. Each point along the boundary of the first Brillouin zone shown as the right inset in Fig. 4(a) determines k_x and k_y in the PBCs in Eq. (19), and thus the matrices in Eqs. (15) and (18). The eigen frequencies of this PC can then be obtained from solving the eigen values of Eqs. (15) and (18). The band diagrams of the TE and TM modes are plotted in Fig. 4(a) and (b), respectively. The circular dots are the results obtained using our compact FDFD algorithm with the index average scheme and the solid lines are obtained using the MIT Photonic-Bands (MPB) package [26] based on the PWE method with 128×128 resolution. It can been seen that the circular dots match quite well with the solid lines for both polarizations. For Fig. 4(a) and (b), 40 grid points are utilized in each lattice distance and it takes just a few minutes to obtain these two figures on Pentium IV 2.0 GHz personal computers.

To examine the efficiency of the index average (IA) scheme, we fix the wavenumber vector k at theMpoint in the first Brillouin zone, i.e., k=(π/a,π/a,0), and search the eigen frequencies using various numbers of grid points with and without the IA scheme. Here and in the following figures, the number of grid points refers to the number of total grid points in the computing domain. In Fig. 5(a) and (b), the lines with circles represent the first eigen frequencies for the TE and TM modes, respectively, which are obtained using our FDFD algorithm without the IA scheme. One can see that it requires quite a lot of grid points to achieve numerical convergence. On the contrary, the lines with squares represent the results using the IA scheme and show quite excellent convergence for both TE and TM modes.

 

Fig. 4. Band diagrams for the 2-D PC formed by square-arranged alumina rods with r/a=0.2 and ε=8.9 in the air. (a) TE mode and (b) TM mode.

Download Full Size | PDF

 

Fig. 5. The first eigen frequency versus the number of grid points for (a) the TE and (b) the TM modes as k is at the M point. The lines with circles and rectangles are the results obtained by our FDFD method without and with the index average scheme, respectively.

Download Full Size | PDF

We then calculate other eigen frequencies for the same wavenumber vector k using our compact FDFD algorithm and three other eigenvalue algorithms for 2-D PCs, which are derived respectively following the FDFD mode solver [20], the FDM proposed by Lüsse et al. [13], and that by Xu et al. [15], with different numbers of total grid points. The relative error here refers to the difference between the result and a reference result obtained with the grid points being 2500. Figure 6(a) and (b) shows the relative errors of the first and second bands for the TE mode versus the number of grid points. All results from the above-mentioned algorithms are calculated with the IA scheme. It shows that as the number of total grid points is larger than 500, all the four algorithms can achieve a relative error less than 1%. Compared with the other three algorithms, our method can provide better convergence behavior thanks to the aid of the compact scheme. As for the TM mode, the relative errors of the first and second bands are presented in Fig. 6(c) and (d), respectively. One can again see the convergence property of our method is good.

 

Fig. 6. 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

As for the 2-D PC with triangular lattice, the band diagrams for both polarizations can be obtained similarly with our method using the modified unit cell. Figure 7(a) and (b) shows the band diagrams for the TE and TM modes of the structure consisting of dielectric cylinders with r=0.2a and relative permittivity ε=11.4 in the air. The triangles are our results using the compact FDFD method with the modified unit cell and the IA scheme, and the solid lines are obtained by the MPB package based on the PWE method [26]. The circles represent the results from the FDTD method based on non-orthogonal coordinate [8] which was proposed to deal with such 2-D PC with non-rectangular unit cell and curved dielectric interfaces. One can see that there exists a visible difference between our results and those from the FDTD method, but fairly good agreement is observed between ours and those by the PWE method. To examine the convergence characteristics, we consider the k vector corresponding to the K point in the first Brillouin zone shown as the right inset in Fig. 7(a), i.e., k=(2π/√3a,2π/3a,0). The convergent behaviors for the first and second bands of the TE and TM modes are illustrated in Fig. 8(a)–(d). Compared with the other three methods as in Fig. 6, our method can still provide better numerical convergence.

We discuss the order of numerical accuracy in the following. Although the compact FD scheme in space is a fourth-order one, as mentioned in [23], an arithmetic average of permittivity at the dielectric interface could reduce the accuracy to second order. In our case, the accuracy of the proposed method is better than second order in the results of Fig. 6(a), (c), and (d) and Fig. 8(a), (c), and (d), whereas the accuracy is a little less than second order in the results of Fig. 6(b) and 8(b). Meade et al. [27] discussed numerical convergence in the PWE calculation of the first band at the zone edge for the square lattice of dielectric rods under different IA schemes. They found that the same IA scheme as in this paper showed very good convergence for the TM band, but not for the TE band. When using an inverse IA scheme involving a weighted average of the inverses of the permittivities, i.e., 1/ε_z(i, j)=f/ε ₁+(1-f)/ε ₂, the convergence performance became reversed between the TM and TE bands. They thus proposed a tensor dielectric scheme which performed well for both bands. In our FDFD calculation, the IA scheme performs well for both first bands as seen in Figs. 6 and 8. We have tried the inverse IA scheme in our FDFD analysis of the TE first and second bands. It reveals that the convergence becomes very good for the TE second band but much worse for the TE first band. Apparently, the effects of the IA schemes on the PWE method and on the FDFD method are different, which is worth further investigation.

 

Fig. 7. Band diagrams of (a) the TE mode and (b) the TM mode for the 2-D PC formed by triangular-arranged dielectric cylinders with r/a=0.2 and ε=11.4 in the air. Our results (triangles) are compared with the results from the FDTD method (circles) and the PWE method (solid lines).

Download Full Size | PDF

 

Fig. 8. The convergence properties of our method for the 2-D triangular-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

4. Conclusion

We have demonstrated a new eigenvalue algorithm for solving the band diagrams of the 2-D PCs, which is based on the finite-difference frequency-domain method associated with the fourth-order accurate compact scheme [23]. The index average scheme is also utilized to deal with the curved dielectric interfaces inside the unit cell. For the non-rectangular unit cell in triangular-lattice PCs, a modified unit cell of rectangular shape is introduced for easier numerical treatment. The band diagrams for both TE and TM modes in 2-D PCs are successfully obtained and good agreement with the PWE method calculation has been achieved. Our method is found to have better convergence property compared with other FD and FDFD algorithms.