Fourier factorization with complex polarization bases in the plane-wave expansion method applied to two-dimensional photonic crystals

We demonstrate an enhancement of the plane wave expansion method treating two-dimensional photonic crystals by applying Fourier factorization with generally elliptic polarization bases. By studying three examples of periodically arranged cylindrical elements, we compare our approach to the classical Ho method in which the permittivity function is simply expanded without changing coordinates, and to the normal vector method using a normal-tangential polarization transform. The compared calculations clearly show that our approach yields the best convergence properties owing to the complete continuity of our distribution of polarization bases. The presented methodology enables us to study more general systems such as periodic elements with an arbitrary cross-section or devices such as photonic crystal waveguides.


Introduction
Photonic crystals (PhCs) are modern artificially structured systems of a broad interest from many viewpoints of fundamental research and applications [1,2,3,4]. Studying their electromagnetic properties is significantly important for developing PhC-based promising applications such as purely optical integrated circuits [5,6,7], artificial metamaterials with high tunability [8,9,10,11], high-sensitivity PhC biosensors [12,13], or devices based on phenomena not accessible in conventional media [14,15,16]. It is also necessary for accounting for structural colors of wings of butterflies or beetles, feathers of birds, or iridescent plants [17,18].
Owing to the PhC periodicity, the plane-wave expansion (PWE) method is conventional for calculating PhC modes and photonic band structures. Its essence is the Fourier expansion of electromagnetic fields and material parameters, which is not only a counterpart of the PWE method for electronic crystals, but it also advantageously follows the classical coupled-wave theory developed for diffraction gratings during last several decades [19,20,21]. One of the most important discoveries of the grating theory are the Fourier factorization (FF) rules for expanding the ratios of the discontinuous permittivity and the discontinuous electric field [22]. The three FF theorems were first formulated for one-dimensional (1D) gratings and then generalized to two-dimensional (2D) systems [23], arbitrary periodic reliefs [24], anisotropic [25] and slanted [26] periodic structures, their various combinations [27,28,29], and other applications [30,31,32]. As will be demonstrated in this paper, the FF rules establish suitable principles for an accurate solution of Maxwell's equations in the plane-wave basis, which considerably enhances the numerical performance compared to previously used implementations.
In the case of 2D-periodic structures, the FF rules were first applied to "zigzag" Fourier expansions [23], which yielded an improvement for rectangular dots or holes, but were not sufficient for the staircase approximation of round elements. After that a coordinate transform was applied to treat individually the normal and tangential components of the electric field on 1D sinusoidal-relief gratings, which enabled the application of the correct rule for each field component [24]. Later David et al. [33] utilized the normal-tangential field separation to 2D PhCs composed of circular and elliptical holes. Similarly, Schuster and colleagues [34] applied this method to 2D gratings, and also suggested more general distributions of polarization bases [35]. These approaches, always dealing with linear polarizations, enabled a significant improvement of the convergence properties, but ignored the fact that the transformation matrix between the Cartesian and the normal-tangential component bases of polarization became discontinuous at the center and along the boundaries of the periodic cell, which slowed down the resulting convergence. To overcome these discontinuities, a distribution of more complex (i.e., generally elliptic) polarization bases was recently suggested to improve optical simulations of 2D gratings [36].
Here we are dealing with 2D PhCs made as 2D-periodic elements of the circular crosssection, for which we implement the PWE method by using FF with complex polarization bases. From the mathematical point of view this complexity only requires complex-valued Jones vectors, so that the generalization from the previous method of David is surprisingly simple. Both methods are compared with respect to their convergence properties, together with the classical method of Ho [3] (where the electric permittivity is expanded without any coordinate transform), to show that the method presented here yields the best performance.

PWE method for 2D structures
We study a 2D PhC composed of infinite cylinders with a circular cross-section with either square [ Fig. 1(a)] or hexagonal [ Fig. 1(b)] periodicity. While the structure is invariant along the z-axis, the periodicity is assumed in the x-and y-axes directions. For the square symmetry the unit cell has the dimensions a x = a y = a, and for the hexagonal symmetry it can be chosen as an Fig. 1. Two studied configurations of 2D PhCs, with square (a) and hexagonal (b) symmetry, with denoted square and rectangular unit cells, respectively. The corresponding first Brillouin zones of the reciprocal space are depicted in (c) and (d), respectively, where the emphasized triangles denote irreducible fractions. Note that for calculations we use a rectangular unit supercell (solid rectangle) rather than a hexagonal primitive cell (dashed hexagon) in the case of the hexagonal lattice (b) so that the corresponding first Brillouin zone is the solid rectangle instead of the dashed hexagon in (d), which creates a folded band structure, from where we choose values corresponding to the usual convention.
a-by-a √ 3 rectangle. The corresponding first Brillouin zones of the reciprocal space are depicted in Figs. 1(c) and 1(d), together with the symmetry points Γ, X, M, and Γ, M, K, respectively. According to this 2D-periodic structure we define and expand the relative permittivity function where ε mn are its Fourier coefficients and G x = 2π/a x and G y = 2π/a y are the reciprocal lattice vectors.
Maxwell's equations for the H-polarization (where the electric field E is in the x-y plane whereas the magnetic field H is along the z-axis) are where ω is the photon frequency, and ε 0 and µ 0 are the permittivity and permeability of vacuum, or Note that here we do not study the E-polarization (for which the electric field E is along the z-axis) because in that case the Ho method satisfies the correct Fourier factorization rules.
Assuming an in-plane anisotropy of the relative permittivity function, defining a scaled electrical displacementD, defining a 2-by-2 matrix of electrical impermittivity η η η = ε ε ε −1 , and substituting Eq. (6) into Eq. (5) yields where η jk are the components of the electrical impermittivity, c = (ε 0 µ 0 ) −1/2 is the light speed in vacuum, and the factor ω 2 /c 2 is calculated as the eigenvalues of the operator on the left-hand side of Eq. (8).
According to the Floquet-Bloch theorem, in the reciprocal space the magnetic field undergoes the Fourier expansion where

Elementary (Cartesian) method (Model A)
In this article we compare several models corresponding to different FF approaches. First, Model A assumes the solution in the basis of the x and y polarizations uniform within the periodic cell [ Fig. 2(a)], where in accordance with the Ho method [3] we choose the Laurent rules The components of the electric impermittivity in Eq. (8) then becomes [[ε]] −1 for the cases of η xx , η yy , and zero for the cases of η xy , η yx , or In Fig. 2(a) we depict the uniform polarization basis within the first quadrant of the square periodic cell, 0 < x < D(0), 0 < y < D(π/2), where D(φ ) is the distance between the cell's center and its boundary (a function of the polar coordinate φ ). We also denote R the radius of the cylindrical element.

Normal vector method (Model B)
According to Li [22] the Laurent rule is valid (for the reason of uniform convergence) for multiplying functions that possess no concurrent discontinuities, whereas the inverse rule is used for functions whose product is continuous. Obviously, neither the Laurent rule nor the inverse rule is correct for both products in Eqs. (10) and (11), because both pairs of functions have concurrent discontinuities and both productsD x andD y are discontinuous as well. On the other hand, by an appropriate change of the polarization bases at all points (using a spacedependent Jones matrix transform F), we can treat independently the normal (u) and tangential (v) components of the fields by the correct rules, The field components E u ,D u are normal to the discontinuities of the relative permittivity function, while E v ,D v are tangential. The FF rules used in Eqs. (14) and (15) where the polar angle φ (x, y) is in the first cell distributed according to the polar coordinates re iφ = x + iy, and then periodically repeated over the entire 2D space. This enables defining the matrices [[c]], [[s]] from the corresponding 2D-periodic functions c = cos φ , s = sin φ . Let u and v be the two columns of the matrix F, both being mutually orthogonal basis vectors of linear polarization [37]. From the above definitions we see that u is a polarization vector normal to the structure discontinuities, whereas v is tangential. In Fig. 2(b) we depict the distribution of the u-v basis within the first quadrant of the periodic cell; the polarization vectors are constant along the lines of the constant azimuth (φ = const) and rotate as φ increases. The distribution of the rotation of u within the square periodic cell is displayed in Fig. 3(a), from where it is obvious that the matrix function F(x, y) has no discontinuities concurrent with the electric field, so that we can use both Laurent and inverse rules for the transformation of polarization, e.g., Combining Eqs. (14), (15), and (17) yields whose components are immediately applicable to Eq. (8).

Method with elliptical polarization bases (Model C)
The above approach (Model B) only deals with linear polarizations and thus suffers from the fact that the matrix function F(x, y) has a singularity at the center of the periodic cell and other discontinuities along the cell boundaries. This slows down the convergence of the numerical implementation, as will be evidenced below.
As suggested in the previous work for the case of diffraction by 2D gratings [36], we can make F continuous by using complex functions ξ and ζ or, in other words, by defining u, v as complex vectors corresponding to generally elliptic polarizations, (which are still orthogonal). By means of rotation θ and ellipticity E we define the first basis vector where Here R denotes the radius of the circular element and is the distance from the cell's center to its edge. In Eq. (22) the Jones vector on the right represents a polarization ellipse (with ellipticity E ) oriented along the x coordinate, the matrix in the middle rotates this polarization by the azimuth θ , and the factor e iθ preserves the continuity of the phase at the center and along the boundaries of the cell. This continuity can be easily checked by evaluating the limits which is the vector of left circular polarization (independent of φ ). The distribution of the u-v polarization basis within the first quadrant of the periodic cell is schematically depicted in Fig. 2(c), where the azimuth of the polarization ellipse is constant along the lines coming from the cell's center, which is similar to Model B; hence, Fig. 3(a) serves for both Model B and Model C. However, the ellipticity is now zero (corresponding to linear polarization) only on the boundaries of the circular element, has the maximum value (π/4 for circular polarization) at the cell's center and along its boundaries, and continuously varies (with a smooth sine dependence) in the intermediate ranges [ Fig. 3(b)]. Finally we obtain a smooth and completely continuous matrix function F(x, y), which is analogously used to calculate the impermittivity in the reciprocal space In the case of the hexagonal periodicity we define u and the other periodic quantities inside one hexagon (half the area of the rectangular unit cell) where we can use formally the same equations as above, except for which is now the distance from the hexagon's center to its edge. Here a is the hexagon's shortest width (equal to the width a x of the rectangular cell). After periodic repeating over the entire 2D space, the azimuth and ellipticity distributions of the basis vector u become those displayed in Figs. 3(e) and 3(f), respectively. Note that although the color distribution in Fig. 3(e) is discontinuous along the hexagon's boundaries, the vector u is actually continuous, because the color difference is only due to the 180 • change. Moreover, the vector u in the limit r → D(φ ) corresponds to left circular polarization, so that its azimuth becomes irrelevant.

Modified method for densely arranged elements (Model C')
To analyze a more complicated situation, we consider a PhC with square periodicity where circular elements are densely arranged near each other, i.e., where the radius R is almost the half width a/2 of the periodic cell. Then the convergence properties of F becomes worse, which affects all the derived quantities. For this reason we again redefine the polarization distribution, as suggested in Fig. 2(d). For the modified Model C' we define u to be still same inside the circle (r < R), but different outside. Assuming the rotation and ellipticity along the boundary of the square cell (where "round" denotes rounding towards the nearest integer), we define the rotation and ellipticity outside the circle (r > R) as Assuming otherwise the same Eqs. (19), (21)

Description of simulations
We examine the numerical performances of all the presented models on three samples of 2D PhCs, for which we calculate the eigenfrequencies ω κ (where the band number κ = 1 stands for the lowest eigenfrequency, κ = 2 for the second lowest, etc.) and the corresponding eigenvectors [H z ] κ of Eq. (8). All convergences will be presented according to the maximum Fourier harmonics retained inside the periodic medium (M, N defined in Appendix), which will be kept same for the x and y directions (M = N). First, Sample S1 is a square array of cylindrical rods of the circular cross-section with the diameter 2R = 500 nm, square period a = 1000 nm, relative permittivity of the rods ε 1 = 9, and relative permittivity of the surrounding medium corresponding to vacuum (ε 2 = 1). For clarity we display the first several eigenmodes at the points of symmetry Γ, X, M in Fig. 4(a). The color scale of the modes corresponds to the scaled amplitude distribution of the magnetic field |H z (x, y)| within one cell of the real space (calculated by Model C with M = N = 12). Each eigenmode is denoted by the symmetry point letter and the band number (e.g., X2 denotes the second band at the point X). By comparing the calculated eigenfrequencies and the amplitude distributions we find that the pairs of eigenmodes Γ3-Γ4 and M2-M3 constitute frequency doublets, i.e., the frequencies for one pair are very close to each other (differences are due to numerical errors). Note that although the amplitudes |H z (x, y)| of the two eigenmodes within a doublet seem identical [e.g., the doublet Γ3-Γ4 in Fig. 4(a)], the two actual eigenmodes, which are complex-valued distributions H z (x, y), are indeed different, being two linearly independent eigenvectors of a degenerate eigenvalue.
Similarly, for Sample S2 we assume exactly the same parameters except the diameter of the rods, now being 2R = 900 nm. This corresponds to densely arranged elements (the distance between two adjacent rods is only 100 nm). The amplitude distributions for the same bands at the same symmetry points are displayed in Fig. 4(b) (calculated by Model C' with M = N = 12), where we again observe the frequency doublets in the same bands.
Finally, for Sample H we consider a hexagonal array of cylindrical holes of the circular crosssection with the diameter 2R = 600 nm, hexagonal periodicity a = 1000 nm (corresponding to the rectangular cell of the dimensions a x = 1 µm, a y = √ 3 µm), relative permittivity of the holes corresponding to vacuum (ε 1 = 1), and relative permittivity of the substrate medium (surrounding holes) ε 2 = 12. Since the lowest (nonzero) eigenfrequencies at the Γ and K points are quite high, we confine ourselves only to the first three eigenmodes at the point M, whose amplitude distributions (within the rectangular cell) are displayed in Fig. 4(c).

Results and discussion
We carefully compared the numerical effectivity of the presented models for all the defined samples, and below we present the convergence results for selected bands.
For Sample S1 the difference among the performances of Models A, B, and C was found as follows: The values of Model B always converge much faster than those of Model A, which is a result expected from previous reports [33]. The values of Model C converge significantly faster than those of Model B for the modes Γ3, Γ4, X1, X4, and M2-M4. This can be explained by the discontinuities of the matrix function F(x, y) along the boundaries and at the center of the periodic cell, which are responsible for the Gibbs phenomenon appearing in the partial Fourier sums in Eq. (20). On the other hand, both Models B and C exhibit similar convergence properties for the modes Γ2, X2, X3, and M1, where the amplitude distribution has almost circular symmetry and is nearly zero near the cell's boundary [ Fig. 4(a)], so that the discontinuities do not manifest themselves.
A few examples of compared convergences are displayed in Fig. 5, where a considerable improvement owing to Model C is visible everywhere except for the Γ2 mode. Nevertheless, according to the scale of the frequency axis at the mode Γ2 no improvement by Model C is necessary since Model B already yields the 5-digits precision. For the modes Γ3, X1, and X4 the precision is improved by at least one order.
For Sample H the convergence of Model A is also much slower than those of Models B and C. In Fig. 6 the values of Model A are even beyond the displayed range of the graphs for the modes M2 and M3. Model C yields the best results for the modes M1 and M2. On the other hand, Models B and C converge similarly for the mode M3; obviously a higher range of N would be required to compare the convergences more adequately.
All the eigenmodes Γ2-Γ4, X1-X4, and M1-M4 of Sample S2 were carefully analyzed with respect to the convergence properties of Models A, B, C, and C'. Here the values of Model C', as expected, converge fastest for all the modes except M1, whereas Model A is again the worst. The eigenmode of the mode M1 is again circularly symmetric and has negligible values near the cell's boundaries [ Fig. 4(b)], so that Model B is not affected by the discontinuities of the polarization transform and yields precise values. On the other hand, Model C is found inappropriate for this sample (yielding even worse convergence than Model B) because the rapid variation of the ellipticity of u near the cell boundaries between two adjacent elements (which are very close to each other) requires more Fourier components than the weak discontinuity of the linear polarization u in Model B. This problem is removed in Model C', as obvious from Figs. 2(d) and 3(d). Three examples of convergences (for the modes X1-X3) are shown in Fig. 7, which clearly justifies the suitability of Model C'. Notice again that the values of Model A are in two cases (modes X1 and X2) beyond the displayed range of the graphs (with error higher by about two orders). In the case of the mode X3 the values of Model A do not differ so much, which can be explained by nearly zero values of the eigenmode X3 at all points of PhC discontinuities, which is obviously not the case of the modes X1 and X2 [ Fig. 4(b)].

Conclusions
We have derived three models with respect to different methods of applying the FF rules and compared their numerical performances. Model A (equal to the Ho method) exhibits the worst convergence properties. Model B gives the best performance in a few cases where the eigenmode is circularly symmetric and has a negligible amplitude near the cell's boundaries (then the problems coupled with polarization discontinuities vanish in zero fields). Model C was found to converge significantly faster than Model B for PhCs without dense arrangement of elements. On the other hand, for dense PhCs, where Model C was found inappropriate, the best results were yielded by Model C', which combines the best properties of Models B and C. We can conclude that the method of FF with complex polarization bases, represented here by Models C and C', produces an enhancement of the PWE method, provided that we choose an appropriate distribution of the polarization basis according to the sample under study.
The method can be straightforwardly generalized to PhCs composed of elements of an arbitrary cross-section. It is possible by replacing the constant value of the element's radius R with a function R(φ ), and by a further generalization of the azimuth θ (r, φ ) and ellipticity E (r, φ ) of the vector u to make it again normal to and linear at the element boundaries (and continuous in the entire 2D space). Moreover, the method can also be used to PhC-based devices such as PhC waveguides by applying the demonstrated methodology to the device supercell (similarly as we were here dealing with the rectangular cell of the hexagonal PhC). It is particularly advantageous for devices where high accuracy is required, e.g., for analyzing defect modes near photonic band edges [38, 39, 40], and for large devices for which the available computer memory enables calculations with only a few Fourier harmonics, e.g., PhC-based fibers with large cladding.