Open Access
Add to BibSonomyAbstract
The moving least-square (MLS) basis is implemented for the real-space band-structure calculation of 2D photonic crystals. A value-periodic MLS shape function is thus proposed in order to represent the periodicity of crystal lattice. Through numerical examples, this MLS method is proved to be a promising scheme for predicting band gaps of photonic crystals.
©2003 Optical Society of America
1. Introduction
The photonic band gap material is nowadays one of central issues in photophysics and optoelectronics communities [1, 2]. Accurate computation of the band structure is thus indispensable for the development of various optoelectronic devices and optical fibers based on these materials [3, 4, 5]. It is well known that the conventional plane wave method has convergence problem arising from the abrupt change in the value of dielectric function across the interface between matrix and inclusion [4, 6]. Special techniques are therefore required for the plane wave method to achieve faster convergence, such as interpolating dielectric function [7] for example. Other numerical methods have also been tried for more efficient calculations of photonic band structures. Among them are the finite-difference time-domain (FDTD) method [8, 9, 10], the multiple-scattering theory [11, 12], the transfer matrix method [13], the finite difference method with effective medium theory [14], the finite element method [15, 16, 17], and most recently, the cell method [18].
In this paper, we apply the novel moving least-square (MLS) method to the band-structure computation of 2D photonic crystals. The MLS method is also known as the meshfree method in the communities of solid and fluid mechanics. It is a real-space scheme, and thus the discontinuity of dielectric function can naturally be modelled without requiring the Fourier transform. In addition, any irregular shape of interface between materials is easily recognized in MLS modelling by simply putting nodes along the interface, which is similar to finite element method. On the other hand, it is a point-based method that does not need grid or mesh in MLS approximation procedure, being different from other real-space techniques. Owing to its distinctive features, the MLS method has widely been applied to various kinds of engineering problems (for the most recent review, refer to [19]). Specifically the method may be well-suited to the adaptive computation, which motivates us to investigate here the feasibility of this method for the band-structure calculation of photonic crystals. While executing the photonic band calculation, the most important issue is how to realize the crystal’s periodicity in the MLS modelling procedure. To this end, we have developed an algorithm introducing a value-periodic MLS shape function. Therefore, the goal of this paper is to demonstrate that this periodic MLS shape function can be applied to highly accurate computations of photonic band structure. We show that its performance is outstanding well beyond the conventional plane wave method throughout numerical examples.
2. The moving least-square (MLS) approximation
The MLS representation of an arbitrary function begins with discretizing a domain by a set of nodes. Nodal value u J is assigned to each node. The approximation u h(x) of a function u(x) at any point x in the domain of computation is then expressed by
where NP is the total number of nodes. By minimizing the localized error residual functional that is expressed by the moving least-square procedure (see Appendix or [19, 20] for details), the MLS shape function NJ (x) is defined as
where p(xJ-x) is the polynomial basis vector and W(xJ-x) the window function. They are here abbreviated as pJ(x)=p(xJ-x) and W J (x)=W(xJ-x). The vector b(x) in Eq. (2) is determined by solving the matrix equation of
where p(0)=[1, 0, …, 0]T and the moment matrix M(x) is defined by
The MLS shape function is highly localized. That is, it vanishes in the most of problem domain, except a small region (so called compact support) of which the center is at the position x. The comprehensive description of MLS-based meshfree methods is found in a recent review paper [19] and references therein.
3. Galerkin formulation and matrix eigen-equations
For the computation of band structures in 2D photonic crystals, we begin with the two scalar equations for solving the unknown field intensity ψ(x) as below,
where λ=(ω/c)2 and ε(x) is the dielectric function. ω is the frequency of the monochromatic electromagnetic wave, c the speed of light, and x the position vector in x-y plane. Applying the Floquet-Bloch theorem, the problems are reduced to solving for u(x) as
where k is the wave vector and u(x) is the function fulfilling the periodicity of lattice structure of which the lattice vector is L, i.e., u(x)=u(x+L).
The Galerkin formulations for Eqs. (7) and (8) result in matrix eigen-equations as,
where, for TE modes,
Fig. 1. The translation-and-searching algorithm for constructing a periodic meshfree shape function. Supports of the shape function are in gray within solid circle. The dashed circle illustrates the rod of photonic crystal.
and, for TM modes,
The upper bars in (10) and (11) imply the complex conjugate, and (x) is the periodic MLS shape function. The periodic nature of shape function is emphasized by the superscript P. We assume for now that the periodic MLS shape function enables the approximated functions to satisfy the periodic boundary conditions required. The derivation of such shape function is given in the next section. The matrix eigen-equations are solved using the eigensystem subroutine package (EISPACK) for the numerical examples in Section 5.
When we discretize the unit cell by a set of nodes (or points), it is geometrically natural to locate nodes exactly on the interface between matrix and inclusion. Unfortunately in this situation, it is not clear which dielectric constant has to be assigned to the interface nodes, because the matrix and inclusion have different dielectric constants from each other. In contrast, the MLS-based meshfree method usually employs a set of integration points, apart from the set of nodes, to calculate the integrals in matrix Eqs. [19]. And the integration points are located in either inclusion or matrix, not on the interface. Therefore, the given value of dielectric constant is straightforwardly assigned to integration points. It makes our implementation very simple, not requiring additional techniques such as, for instance, averaging or smoothing the dielectric function.
4. Value-periodic MLS shape functions
Now we modify the MLS representation in order to make the shape function inherently suitable for periodicity. For the computation of the shape function N J (x) at a specific point x, we have to search the associated nodes denoted by J due to the localized property of the shape function. These nodes are located under the support in which the shape function does not vanish. While searching those nodes, we must take the periodicity into consideration. This can be implemented step by step, as shown in Fig. 1 for the point x near the lower left corner of the parallelogram. We first find nodes (inside the unit cell) associated with the original point x. Secondly, we translate the point x to x+a1 and search interior nodes again. Next, it is repeated for x→x+a2 and for x→x+a1+a2. This translation-and-searching procedure is carried out for all translated points at which the support of shape function covers the unit cell. The resulting supports are illustrated in gray color in Fig. 1.
The description is now generalized to an arbitrary point x in the unit cell. Using the lattice vector L=n 1a1+n 2a2, the MLS approximation can be written in general form as
where N J (x+L) is defined by
Note that b is a function of x, not of x+L. The vector b(x) is obtained by solving
The summation with respect to L is performed over the current unit cell itself and the neighboring cells only, because, for any farther L, the corresponding support does not overlap the unit cell of computation. In conclusion, the 2D value-periodic MLS shape function for a parallelogram unit cell whose primitive vectors are a1 and a2, can be written as
Therefore the summation over n i’s results in 9 cells in 2D (or 27 cells for 3D case) at most. That is, n i=-1, 0, +1 in Eq. (15) is enough for any x in the unit cell. It should be carefully noted that the opposing boundaries are physically identical, i.e., Γ1=Γ3 and Γ2=Γ4 as illustrated in Fig. 1. Therefore, for the nodal summation over J=1, …, NP, we must involve one boundary only among each identical pair in order to avoid over-summing. In other word, the nodes on Γ3 and Γ4 have to make no contribution to the summation with respect to J, once the corresponding nodes along Γ1 and Γ2 are summed up. An example of 2D periodic MLS shape function is plotted in Fig. 2 for a parallelogram unit cell.
5. Numerical examples
For numerical examples in this section, employed were the quartic spline function for window function WJ (x), and the two-dimensional linear basis vector of p(x)=[1, x, y]T when calculating the MLS shape function. In general, window function is localized, nonnegative, continuous, and monotonously decreasing from its center. Localization implies that window function WJ (x) does not vanish only in the compact support around its center. The size of support can be controlled by introducing a dilation parameter ρ in defining window function as WJ (x; ρ). The compactness is necessary for efficiency. Once these requirements are met, it is not significant which actual type of window function is employed. From these points of view, quartic spline is one of the common choices for window function in MLS-based meshfree methods. The selection of basis vector p(x) is closely related to convergence. Equations (5) and (6), the master equations, are second order differential equations. In this case, it has been known that the linear basis vector is the simplest guaranteeing convergence [19]. For photonic band gap materials, effects of higher order basis vectors on the accuracy of band-structure computations may be one of the future studies.
5.1 Electromagnetic Kronig-Penney problem
The electromagnetic Kronig-Penney problem is now considered. The unit cell is depicted in Fig. 3, which is composed of two different dielectric materials. The analytic solution of this case can be obtained elsewhere [17]. This example is taken into consideration in order to investigate how accurately the MLS method handles the discontinuity of dielectric function. Dielectric constants of ε1=1.0 and ε2=9.0 are used for the simulation. Frequency band structures of both TM and TE modes are given in Fig. 4. Lines are for the analytic solution and open circles for the MLS results using 11×11 nodes for TM modes and 41×41 nodes for TE modes. Convergence rates of the lowest five eigenvalues at point of k=(π/a)[0.5, 0.5] are given in Fig. 5 in which En ’s denote the normalized frequency eigenvalues (ωn a/2πc) and Ee is the analytic solution. The average slope of TM and TE modes are 3.12 and 2.77 respectively. The results for TM modes show better performance over the TE modes case, not only in convergence but also in accuracy. This is due to the difference of smoothness properties of eigenvectors [16]. There may be several ways to enhance the accuracy of MLS results on TE modes. The meshfree treatment of discontinuity using Lagrange multiplier has been applied to solid mechanics by [21]. Its application to photonic crystals remains for our future study.
5.2 Square lattice of circular rods
Photonic band structures of square lattice comprised of circular rods are computed in this section. The square lattice structure and the square unit cell are illustrated in Fig. 6. The unit cell in the figure is discretized by 1697 nodes. The radius of circular cross section is r=0.2a for this example. Dielectric constants are εm=1.0 and εr=8.9 for matrix and rods respectively. In this case, there exists a wide band gap in TM modes as widely known in literatures. Our results on band structures also confirm the wide band gap in TM modes as shown in Fig. 7. In this example, the width of band gap is determined by the 1st TM mode at point M and the 2nd TM mode at point X. The numerical results at these points obtained by both the MLS method and plane wave method are compared in Fig. 8, where the faster converging behavior of the method over the plane wave method are demonstrated. In the case of MLS method, the x-axis indicates the size of discretization such as h=a/(NP -1), where NP denotes the number of nodes in one direction. For the plane wave method, it similarly implies h=a/(NPW -1), where NPW means the number of plane wave used for the computation. This notation is arbitrary because the MLS method solves generalized eigenvalue problems, while the plane wave method does standard eigenvalue problems. The MLS method thus requires more computing time than the conventional plane wave method in general. Nonetheless, it is convenient for the purpose of direct comparison as in Fig. 8, because both (NP)2 and (NPW )2 determine the size of system matrix to be solved.
Fig. 4. Band structures of the electromagnetic Kronig-Penney problem. MLS results are denoted by open circles and analytic results by solid lines.
Fig. 5. Convergence rates of the lowest five eigenvalues for the electromagnetic Kronig-Penney problem at k=(π/a)(0.5, 0.5).
Fig. 8. Comparisons of MLS method and plane wave method (TM modes) for the square lattice of circular rods.
6. Concluding remarks
In this paper, the value-periodicMLS shape function is implemented in order to calculate the frequency band-structure of 2D photonic crystals. Through numerical examples, it has been verified that MLS results are more accurate than those of the conventional plane wave methods. Specifically TM modes results of meshfree method are outstanding, compared with TE mode results. This property has also been found in other real-space techniques such as finite element method and finite difference method [14, 16]. In order to enhance the accuracy and convergence of TE modes solutions by the method, special techniques may be applied in our future study. Among them is the meshfree treatment of material discontinuity using Lagrange multipliers, which has been developed for solid mechanics by [21]. Extending MLS band-structure computation to 3D photonic crystal is an example to which our periodic MLS shape function can immediately be applied. The effect of mechanical deformation on the change of photonic band is also an interesting subject.
Appendix
Here we describe a typical procedure how to define the shape function of the moving least-square (MLS) method. Although there are other routes to defining it, we mainly follow a recent Ref. [20]. The MLS method begins with the local approximation of an arbitrary function u(x) as
in a small region around x̄. The column vector p(x) must contain at least 1 and monomials in order to fulfill the linear reproducing conditions. For example, p(x)=[1, x, y, z]T for 3D. By reproducing conditions, we mean that the MLS approximation can reproduce exactly the functions which appears in p(x). This is the well-known feature of MLS-based meshfree methods. Any polynomials of higher order or arbitrary types of functions can be added in the component of p(x) for the specific purpose of enrichment in MLS approximation.
After discretizing a problem domain by a set of points (I=1, …, NP), the idea of moving least-square interpolant is used in determining the coefficient vector a(x̄). That is, a(¯x) is determined by minimizing the local error residual functional that is expressed as
where the window function W(x) is a compactly-supported continuous function. By solving ∂J/∂a=0 for a(x̄), we have
where the moment matrix M(x̄) is defined by
Next we insert Eq. (18) into Eq. (16) for a(x̄), and then obtain the best local approximation of u(x) as
Finally, the global MLS approximation of u(x) is obtained by taking the limit of Eq. (20) as x̄→x,
where p(0)=[1, 0, …, 0]T. The moment matrix M(x) becomes
Equation (21) can now be rewritten by given sampling values u I=u(xI) and the MLS shape function N I (x)=N(xI-x) as
where
Due to the symmetry of M(x), it is straightforward to verify that the shape function can also be written as
which is identical to Eq. (2) in Section 2. For comprehensive and systematic derivations of MLS shape function and its derivatives, it is recommended to refer to Refs. [19, 20].
Acknowledgements
This research was supported by a grant (02-K14-01-012-1-0) from Center for Nanoscale Mechatronics & Manufacturing of 21st Century Frontier Research Program of Korea.
References and links
1. J.D. Joannopoulos, R.D. Meade, and J.N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton University Press, Princeton, 1995).
2. Y. Xia, “Photonic crystals,” Adv. Mater. 13, 369 (2001) and papers in this special issue. [CrossRef]
3. K. Busch, “Photonic band structure theory: assesment and perspectives,” C. R. Physique 3, 53–66 (2002). [CrossRef]
4. D. Cassagne, “Photonic band gap materials,” Ann. Phys. Fr. 23(4), 1–91 (1998). [CrossRef]
5. J.B. Pendry, “Calculating photonic band structure,” J. Phys.: Condens. Matter 8, 1085–1108 (1996). [CrossRef]
6. H.S. Sözüer, J.W. Haus, and R. Inguva, “Photonic bands: Convergence problems with the plane-wave method,” Phys. Rev. B 45, 13962–13972 (1992). [CrossRef]
7. R.D. Meade, A.M. Rappe, K.D. Brommer, J.D. Joannopoulos, and O.L. Alerhand, “Accurate theoretical analysis of photonic band-gap materials,” Phys. Rev. B 48, 8434–8437 (1993). [CrossRef]
8. 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]
9. A.J. Ward and J.B. Pendry, “Calculating photonic Green’s functions using a nonorthogonal finite-difference time-domain method,” Phys. Rev. B 58, 7252–7259 (1998). [CrossRef]
10. 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,” J. Appl. Phys. 87, 8268–8275 (2000). [CrossRef]
11. K.M. Leung and Y. Qiu, “Multiple-scattering calculation of the two-dimensional photonic band structure,” Phys. Rev. B 48, 7767–7771 (1993). [CrossRef]
12. X. Wang, X.G. Zhang, Q. Yu, and B.N. Harmon, “Multiple-scattering theory for electromagnetic waves,” Phys. Rev. B 47, 4161–4167 (1993). [CrossRef]
13. J.B. Pendry and A. MacKinnon, “Calculation of photon dispersion relations,” Phys. Rev. Lett. 69, 2772–2775 (1992). [CrossRef] [PubMed]
14. L. Shen, S. He, and S. Xiao, “A finite-diference eigenvalue algorithm for calculating the band structure of a photonic crystal,” Comput. Phys. Comm. 143, 213–221 (2002). [CrossRef]
15. W. Axmann and P. Kuchment, “An efficient finite element method for computing spectra of photonic and acoustic band-gap materials: I. Scalar case,” J. Comput. Phys. 150, 468–481 (1999). [CrossRef]
16. D.C. Dobson, “An efficient band structure calculations in 2D photonic crystals,” J. Comput. Phys. 149, 363–376 (1999). [CrossRef]
17. C. Mias, J.P. Webb, and R.L. Ferrari, “Finite element modelling of electromagnetic waves in doubly and triply periodic structures,” IEE Proc.-Optoelectron. , 146(2), 111–118 (1999). [CrossRef]
18. M. Marrone, V.F. Rodriguez-Esquerre, and H.E. Hernandez-Figueroa, “Novel numerical method for the analysis of 2D photonic crystals: the cell method,” Opt. Express 10, 1299–1304 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-22-1299 [CrossRef] [PubMed]
19. S. Li and W.K. Liu, “Meshfree and particle methods and their applications,” Applied Mechanics Review , 55, 1–34 (2002). [CrossRef]
20. D.W. Kim and Y. Kim, “Point collocation method using the fast moving least-square reproducing kernel approximation,” Int. J. Numer. Methods Engrg. 56, 1445–1464 (2003). [CrossRef]
21. L.W. Cordes and B. Moran, “Treatment of material discontinuity in the Element-Free Galerkin method,” Comput. Meth. Appl. Mech. Eng. 139, 75–89 (1996). [CrossRef]
