Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

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

Open Access Open Access

Abstract

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.

© 2010 Optical Society of America

1. Introduction

Photonic crystals (PhCs) are modern artificially structured systems of a broad interest from many viewpoints of fundamental research and applications [14]. Studying their electromagnetic properties is significantly important for developing PhC-based promising applications such as purely optical integrated circuits [57], artificial metamaterials with high tunability [811], high-sensitivity PhC biosensors [12, 13], or devices based on phenomena not accessible in conventional media [1416]. 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 [1921]. 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 [2729], and other applications [3032]. 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 cross-section, 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.

2. 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 ax = ay = a, and for the hexagonal symmetry it can be chosen as an a-by-a3 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

ɛ(x,y)=m,n=+ɛmnei(mGxx+nGyy),
where ɛmn are its Fourier coefficients and Gx = 2π/ax and Gy = 2π/ay are the reciprocal lattice vectors.

 figure: Fig. 1

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.

Download Full Size | PDF

Maxwell’s equations for the H-polarization (where the electric field E is in the xy plane whereas the magnetic field H is along the z-axis) are

yHz=iωɛ0ɛEx,
xHz=iωɛ0ɛEy,
xEyyEx=iωμ0Hz,
where ω is the photon frequency, and ɛ0 and μ0 are the permittivity and permeability of vacuum, or
[y,x][ExEy]=iωμ0Hz,
ɛ[ExEy]=1iωɛ0[yx]Hz.
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 displacement ,

[D˜xD˜y]=ɛ[ExEy]=[ɛxxɛxyɛyxɛyy][ExEy],
defining a 2-by-2 matrix of electrical impermittivity η = ɛ−1, and substituting Eq. (6) into Eq. (5) yields
(yηxxy+xηyxy+yηxyxxηyyx)Hz=ω2c2Hz,
where ηjk are the components of the electrical impermittivity, c = (ɛ0μ0)−1/2 is the light speed in vacuum, and the factor ω2/c2 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

Hz(x,y)=m,n=+hz,mnei(pmx+qny),
where pm = kx + mGx, qn = ky + nGy. Analogously, the Fourier components of Ex, Ey, x, y are denoted ex,mn, ey,mn, x,mn, y,mn. Then, as demonstrated in Appendix, the operators of the partial derivatives x, y become the diagonal matrices −ip, −iq, the operators of multiplication by functions such as ɛxx or ηxx become matrices denoted [[ɛxx]] or [[ηxx]], and the field components Ex, Ey, Hz become column vectors denoted [Ex], [Ey], [Hz]. Using these definitions, the multiplications in Eq. (7) can be treated either by the Laurent FF rule, [ɛxxEx] = [[ɛxx]][Ex] etc., or by the inverse FF rule, [ɛxxEx] = [[1/ɛxx]]−1[Ex.]

3. Methods of Fourier factorization

3.1. 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

[D˜x]=[ɛEx]=[[ɛ]][Ex],
[D˜y]=[ɛEy]=[[ɛ]][Ey].
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
[[η]]A=[[[ɛ]]1[[0]][[0]][[ɛ]]1].
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.

 figure: Fig. 2

Fig. 2 Schematic description of the polarization distributions in Models A (a), B (b), C (c), and C’ (d) within the first quadrant of the periodic cell, as functions of the polar coordinates re = x + iy. In (a) the xy Cartesian basis is uniform; in (b) the polarization vectors are normal (u) and tangential (v) to the cylindrical element and constant along lines coming from the center; in (c) and (d) the u, v polarizations are in general elliptic which enables their continuity.

Download Full Size | PDF

3.2. 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 products x and y are discontinuous as well. On the other hand, by an appropriate change of the polarization bases at all points (using a space-dependent Jones matrix transform F),

[ExEy]=F[EuEv],
we can treat independently the normal (u) and tangential (v) components of the fields by the correct rules,
[D˜u]=[[1/ɛ]]1[Eu],
[D˜v]=[[ɛ]][Ev].
The field components Eu, u are normal to the discontinuities of the relative permittivity function, while Ev, v are tangential. The FF rules used in Eqs. (14) and (15) are justified simply because Ev and u are continuous.

A suitable distribution of the matrix F within the periodic cell can obviously be the rotation [33]

F=[cosϕsinϕsinϕcosϕ],
where the polar angle ϕ (x,y) is in the first cell distributed according to the polar coordinates re = 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 uv 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.,

[[Ex][Ey]]=[[F]][[Eu][Ev]],
[[F]]=[[[c]][[s]][[s]][[c]]].
Combining Eqs. (14), (15), and (17) yields
[[Ex][Ey]]=[[F]][[[1ɛ]][[0]][[0]][[ɛ]]1][[F1]][[D˜x][D˜y]],
from where we derive the electric impermittivity in the reciprocal space (corresponding to Model B)
[[η]]B=[[F]][[[1ɛ]][[0]][[0]][[ɛ]]1][[F1]]=[[[1ɛ]][[c2]]+[[ɛ]]1[[s2]],[[1ɛ]][[cs]][[ɛ]]1[[cs]][[1ɛ]][[cs]][[ɛ]]1[[cs]],[[1ɛ]][[s2]]+[[ɛ]]1[[c2]]],
whose components are immediately applicable to Eq. (8).

 figure: Fig. 3

Fig. 3 Distribution of the rotation and ellipticity of the basis polarization vector u for the presented models. The rotation for Models B and C is in (a) [(e) for the hexagonal periodicity]; the ellipticity for Model C is in (b) [(f) for the hexagonal periodicity]; and the rotation and ellipticity for Model C’ are in (c) and (d), respectively. The color scales for both the rotation and the ellipticity are on the right (in degrees). The structure discontinuities (the circular boundaries of the periodic elements) are plotted as black or white circles. Notice that u is always linear along and normal to the circle.

Download Full Size | PDF

3.3. 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,

u=[ξζ],v=[ζ*ξ*]
(which are still orthogonal). By means of rotation θ and ellipticity ℰ we define the first basis vector
u=eiθ[cosθsinθsinθcosθ][cosisin],
where
θ(r,ϕ)=ϕ,
(r,ϕ)={π8(1+cosπrR)(rR)π8{1+cosπ[r+D(ϕ)2R]D(ϕ)R}(r>R).
Here R denotes the radius of the circular element and
D(ϕ)=a/2max(|cosϕ|,|sinϕ|)
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 ℰ) oriented along the x coordinate, the matrix in the middle rotates this polarization by the azimuth θ, and the factor eiθ 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
limr0u=limrD(ϕ)u=12[1i],
which is the vector of left circular polarization (independent of ϕ).

The distribution of the uv 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

[[η]]C=[[[1ɛ]][[ξξ*]]+[[ɛ]]1[[ζζ*]],[[1ɛ]][[ξζ*]][[ɛ]]1[[ξζ*]][[1ɛ]][[ξ*ζ]][[ɛ]]1[[ξ*ζ]],[[1ɛ]][[ζζ*]]+[[ɛ]]1[[ξξ*]]].

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

D(ϕ)=a/2maxn=0,,5[cos(ϕnπ3)],
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 ax 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 rD(ϕ) corresponds to left circular polarization, so that its azimuth becomes irrelevant.

3.4. 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

θb(ϕ)=θ(D(ϕ),ϕ)=π2round(ϕ/π2),
b(ϕ)=(D(ϕ),ϕ)=π8(1cos4ϕ)
(where “round” denotes rounding towards the nearest integer), we define the rotation and ellipticity outside the circle (r > R) as
θ(r,ϕ)=12{θb(ϕ)+ϕ+[θb(ϕ)ϕ]cosπ[r+D(ϕ)2R]D(ϕ)R},
(r,ϕ)=b(ϕ)2{1+cosπ[r+D(ϕ)2R]D(ϕ)R}.
Assuming otherwise the same Eqs. (19), (21), (22), and (25), we obtain for [[η]]C′ formally the same matrix as in Eq. (27), except that the functions ξ and ζ are now derived from different azimuth and ellipticity distributions of u [displayed in Figs. 3(c) and 3(d), respectively]. Note that u is again continuous along the cell’s boundaries; to evaluate its precise limits [when x → ±D(0), y = const or y → ±D(ϕ/2), x = const] is now more complicated.

4. Numerical examples

4.1. 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 [Hz]κ 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 |Hz(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 |Hz(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 Hz(x, y), are indeed different, being two linearly independent eigenvectors of a degenerate eigenvalue.

 figure: Fig. 4

Fig. 4 Amplitude distribution of the scaled magnetic field |Hz(x, y)| within one cell of the real space for Samples S1 (a), S2 (b) and H (c). The color scale, same for each sample, is in the top left corner of (a). The subfigures represent eigenmodes distinguished by the symmetry point letters (Γ, X, M) and band numbers (1–4). The structure discontinuities (boundaries of the rods or holes) are plotted as white or black circles.

Download Full Size | PDF

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 cross-section with the diameter 2R = 600 nm, hexagonal periodicity a = 1000 nm (corresponding to the rectangular cell of the dimensions ax = 1 μm, ay=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).

4.2. 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.

 figure: Fig. 5

Fig. 5 Convergence properties of normalized eigenfrequencies (ωa/2πc) calculated for selected bands of Sample S1. The modes Γ2, Γ3 are on the left part of the figure; X1, X4 in the middle; and M2, M4 on the right. The Models A, B, and C are compared to each other, plotted as crosses, squares, and circles, respectively.

Download Full Size | PDF

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.

 figure: Fig. 6

Fig. 6 Convergence properties of normalized eigenfrequencies (ωa/2πc) calculated for the modes M1 (left), M2 (middle), and M3 (right) of Sample H. The Models A, B, and C are plotted as crosses, squares, and circles, respectively.

Download Full Size | PDF

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)].

 figure: Fig. 7

Fig. 7 Convergence properties of normalized eigenfrequencies (ωa/2πc) calculated for the modes X1 (left), X2 (middle), and X3 (right) of Sample S2. The Models A, B, C, and C’ are plotted as crosses, squares, circles, and triangles, respectively.

Download Full Size | PDF

5. 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 eigen-mode 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 ℰ(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 [3840], 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.

Appendix: Matrices for 2D structures

Here we follow the matrix formation procedures and notation as explained by Visnovsky and Yasumoto [41]. Let f(x,y) be one component of the pseudo-periodic electric field inside the medium described by the 2D-periodic function ɛ(x,y). Their Fourier expansions are written

ɛ(x,y)=m,n=+ɛmnei(mGxx+nGyy),
f(x,y)=m,n=+fmnei(pmx+qny).
We will here derive the matrix expressions of the fundamental relations
h(x,y)=ɛ(x,y)f(x,y),
gx(x,y)=xf(x,y),
gy(x,y)=yf(x,y),
i.e., the relations of multiplication by a function and applying partial derivatives. Assuming the expansions of the new functions
h(x,y)=m,n=+hmnei(pmx+qny),
gx(x,y)=m,n=+gx,mnei(pmx+qny),
gy(x,y)=m,n=+gy,mnei(pmx+qny),
we rewrite Eqs. (35)(37) using the convolution rule, and applying the partial derivatives as follows:
hmn=k,l=+ɛmk,nlfkl,
gx,mn=ipmfmn,
gy,mn=iqnfmn.
Assuming furthermore a finite number of the retained Fourier coefficients, i.e., using the summation m=M+Mn=N+N, we can renumber all the indices so that instead of a couple of two sets m ∈ {−M, −M + 1, ..., M} and n ∈ {−N, −N + 1, ..., N} we have a single set of indices α ∈ {1, 2, ..., αmax}, with αmax = (2M + 1)(2N + 1), related
α(m,n)=m+M+1+(n+N)(2M+1),
n(α)=(α1)div(2M+1)N,
m(α)=(α1)mod(2M+1)M,
where “div” denotes the operation of integer division and “mod” the remainder (the modulo operation). Then we can rewrite Eqs. (41)(43) into the matrix relations
[h]=[[ɛ]][f],
[gx]=ip[f],
[gy]=iq[f],
where [f], [h], [gx,] and [gy] are column vectors whose αth elements are the Fourier [m,n] elements of the functions f, h, gx, and gy, indexed by α (m, n) defined in Eq. (44), and where [[ɛ]], p, and q are matrices whose elements are defined
[[ɛ]]αβ=ɛm(α)m(β),n(α)n(β),
pαβ=pm(α)δαβ,
qαβ=qn(α)δαβ,
where the indices on the right hand parts are defined by Eqs. (45) and (46) and where δαβ denotes the Kronecker delta. As a summary we can say that the multiplication by a function is in the reciprocal space represented by the matrix [[ɛ]] (in the sense of the limit αmax → ∞), which is a generalization of the Toeplitz matrix used for 1D periodicity, and that the partial derivatives are represented by the diagonal matrices −ip and −iq.

Acknowledgments

This work is part of the research plan MSM 0021620834 financed by the Ministry of Education of the Czech Republic and was partially supported by a Marie Curie International Reintegration Grant (no. 224944) within the 7th European Community Framework Programme and by the Grant Agency of the Czech Republic (no. 202/09/P355 and P204/10/P346). One of the authors (R.A.) thanks Koki Watanabe and Mathias Vanwolleghem for fruitful discussions.

References and links

1. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton Univ., 1995).

2. E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. 58, 2059–2062 (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. J. D. Joannopoulos, P. R. Villeneuve, and S. Fan, “Photonic crystals: Putting a new twist on light,” Nature 386, 143–149 (1997). [CrossRef]  

5. S.-Y. Lin, E. Chow, V. Hietala, P. R. Villeneuve, and J. D. Joannopoulos, “Experimental demonstration of guiding and bending of electromagnetic waves in a photonic crystal,” Science 282, 274–276 (1998). [CrossRef]   [PubMed]  

6. S. Noda, “Recent progresses and future prospects of two- and three-dimensional photonic crystals,” J. Lightwave Technol. 24, 4554–4567 (2006). [CrossRef]  

7. J. P. Hugonin, P. Lalanne, T. P. White, and T. F. Krauss, “Coupling into slow-mode photonic crystal waveguides,” Opt. Lett. 32, 2638–2640 (2007). [CrossRef]   [PubMed]  

8. S. Datta, C. T. Chan, K. M. Ho, and C. M. Soukoulis, “Effective dielectric constant of periodic composite structures,” Phys. Rev. B 48, 14936–14943 (1993). [CrossRef]  

9. F. Genereux, S. W. Leonard, H. M. van Driel, A. Birner, and U. Gosele, “Large birefringence in two-dimensional silicon photonic crystals,” Phys. Rev. B 63, 161101 (2001). [CrossRef]  

10. A. A. Krokhin, P. Halevi, and J. Arriaga, “Long-wavelength limit (homogenization) for two-dimensional photonic crystals,” Phys. Rev. B 65, 115208 (2002). [CrossRef]  

11. E. Reyes, A. A. Krokhin, and J. Roberts, “Effective dielectric constants of photonic crystal of aligned anisotropic cylinders and the optical response of a periodic array of carbon nanotubes,” Phys. Rev. B 72, 155118 (2005). [CrossRef]  

12. N. Skivesen, A. Tetu, M. Kristensen, J. Kjems, L. H. Frandsen, and P. I. Borel, “Photonic-crystal waveguide biosensor,” Opt. Express15, 3169–3176 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-6-3169. [CrossRef]   [PubMed]  

13. I. D. Block, N. Ganesh, M. Lu, and B. T. Cunningham, “Sensitivity model for predicting photonic crystal biosensor performance,” IEEE Sens. J. 8, 274–280 (2008). [CrossRef]  

14. H. Kosaka, T. Kawashima, A. Tomita, M. Notomi, T. Tamamura, T. Sato, and S. Kawakami, “Superprism phenomena in photonic crystals,” Phys. Rev. B 58, R10096–R10099 (1998). [CrossRef]  

15. A. A. Krokhin and E. Reyes, “Homogenization of magnetodielectric photonic crystals,” Phys. Rev. Lett. 93, 023904 (2004). [CrossRef]   [PubMed]  

16. H. Benisty, “Dark modes, slow modes, and coupling in multimode systems,” J. Opt. Soc. Am. B 26, 718–724 (2009). [CrossRef]  

17. P. Vukusic and J. R. Sambles, “Photonic structures in biology,” Nature 424, 852–855 (2003). [CrossRef]   [PubMed]  

18. S. Kinoshita and S. Yoshioka, “Structural colors in nature: The role of regularity and irregularity in the structure,” ChemPhysChem 6, 1442–1459 (2005). [CrossRef]   [PubMed]  

19. R. Petit (ed.), Electromagnetic Theory of Gratings (Springer, Berlin, 1980).

20. M. Neviere and E. Popov, Light Propagation in Periodic Media: Diffraction Theory and Design (Marcel Dekker, New York, 2003).

21. D. Maystre, “Rigorous vector theories of diffraction gratings,” Prog. Opt. 21, 1–67 (1984). [CrossRef]  

22. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13, 1870–1876 (1996). [CrossRef]  

23. L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14, 2758–2767 (1997). [CrossRef]  

24. E. Popov and M. Neviere, “Grating theory: new equations in Fourier space leading to fast converging results for TM polarization,” J. Opt. Soc. Am. A 17, 1773–1784 (2000). [CrossRef]  

25. L. Li, “Reformulation of the Fourier modal method for surface-relief gratings made with anisotropic materials,” J. Mod. Opt. 45, 1313–1334 (1998). [CrossRef]  

26. B. Chernov, M. Neviere, and E. Popov, “Fast Fourier factorization method applied to modal analysis of slanted lamellar diffraction gratings in conical mountings,” Opt. Commun. 194, 289–297 (2001). [CrossRef]  

27. K. Watanabe, R. Petit, and M. Neviere, “Differential theory of gratings made of anisotropic materials,” J. Opt. Soc. Am. A 19, 325–334 (2002). [CrossRef]  

28. K. Watanabe, “Numerical integration schemes used on the differential theory for anisotropic gratings,” J. Opt. Soc. Am. A 19, 2245–2252 (2002). [CrossRef]  

29. L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A 5, 345–355 (2003).

30. P. Boyer, E. Popov, M. Neviere, and G. Tayeb, “Diffraction theory in TM polarization: application of the fast Fourier factorization method to cylindrical devices with arbitrary cross section,” J. Opt. Soc. Am. A 21, 2146–2153 (2004). [CrossRef]  

31. N. Bonod, E. Popov, and M. Neviere, “Light transmission through a subwavelength microstructured aperture: electromagnetic theory and applications,” Opt. Commun. 245, 355–361 (2005). [CrossRef]  

32. N. Bonod, E. Popov, and M. Neviere, “Fourier factorization of nonlinear Maxwell equations in periodic media: application to the optical Kerr effect,” Opt. Commun. 244, 389–398 (2005). [CrossRef]  

33. A. David, H. Benisty, and C. Weisbuch, “Fast factorization rule and plane-wave expansion method for two-dimensional photonic crystals with arbitrary hole-shape,” Phys. Rev. B 73, 075107 (2006). [CrossRef]  

34. T. Schuster, J. Ruoff, N. Kerwien, S. Rafler, and W. Osten, “Normal vector method for convergence improvement using the RCWA for crossed gratings,” J. Opt. Soc. Am. A 24, 2880–2890 (2007). [CrossRef]  

35. P. Gotz, T. Schuster, K. Frenner, S. Rafler, and W. Osten, “Normal vector method for the RCWA with automated vector field generation,” Opt. Express16, 17295–17301 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-22-17295. [CrossRef]   [PubMed]  

36. R. Antos, “Fourier factorization with complex polarization bases in modeling optics of discontinuous bi-periodic structures,” Opt. Express17, 7269–7274 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-9-7269. [CrossRef]   [PubMed]  

37. R. M. A. Azzam and N. M. Bashara, Ellipsometry and Polarized Light (North Holland, Amsterdam, 1997).

38. S. Mahmoodian, C. G. Poulton, K. B. Dossou, R. C. McPhedran, L. C. Botten, and C. M. de Sterke, “Modes of shallow photonic crystal waveguides: semi-analytic treatment,” Opt. Express17, 19629–19643 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-22-19629. [CrossRef]   [PubMed]  

39. S. Mahmoodian, R. C. McPhedran, C. M. de Sterke, K. B. Dossou, C. G. Poulton, and L. C. Botten, “Single and coupled degenerate defect modes in two-dimensional photonic crystal band gaps,” Phys. Rev. A 79, 013814 (2009). [CrossRef]  

40. K. B. Dossou, C. G. Poulton, L. C. Botten, S. Mahmoodian, R. C. McPhedran, and C. M. de Sterke, “Modes of symmetric composite defects in two-dimensional photonic crystals,” Phys. Rev. A 80, 013826 (2009). [CrossRef]  

41. S. Visnovsky and K. Yasumoto, “Multilayer anisotropic bi-periodic diffraction gratings,” Czech. J. Phys. 51, 229–247 (2001). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1
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.
Fig. 2
Fig. 2 Schematic description of the polarization distributions in Models A (a), B (b), C (c), and C’ (d) within the first quadrant of the periodic cell, as functions of the polar coordinates re = x + iy. In (a) the xy Cartesian basis is uniform; in (b) the polarization vectors are normal (u) and tangential (v) to the cylindrical element and constant along lines coming from the center; in (c) and (d) the u, v polarizations are in general elliptic which enables their continuity.
Fig. 3
Fig. 3 Distribution of the rotation and ellipticity of the basis polarization vector u for the presented models. The rotation for Models B and C is in (a) [(e) for the hexagonal periodicity]; the ellipticity for Model C is in (b) [(f) for the hexagonal periodicity]; and the rotation and ellipticity for Model C’ are in (c) and (d), respectively. The color scales for both the rotation and the ellipticity are on the right (in degrees). The structure discontinuities (the circular boundaries of the periodic elements) are plotted as black or white circles. Notice that u is always linear along and normal to the circle.
Fig. 4
Fig. 4 Amplitude distribution of the scaled magnetic field |Hz(x, y)| within one cell of the real space for Samples S1 (a), S2 (b) and H (c). The color scale, same for each sample, is in the top left corner of (a). The subfigures represent eigenmodes distinguished by the symmetry point letters (Γ, X, M) and band numbers (1–4). The structure discontinuities (boundaries of the rods or holes) are plotted as white or black circles.
Fig. 5
Fig. 5 Convergence properties of normalized eigenfrequencies (ωa/2πc) calculated for selected bands of Sample S1. The modes Γ2, Γ3 are on the left part of the figure; X1, X4 in the middle; and M2, M4 on the right. The Models A, B, and C are compared to each other, plotted as crosses, squares, and circles, respectively.
Fig. 6
Fig. 6 Convergence properties of normalized eigenfrequencies (ωa/2πc) calculated for the modes M1 (left), M2 (middle), and M3 (right) of Sample H. The Models A, B, and C are plotted as crosses, squares, and circles, respectively.
Fig. 7
Fig. 7 Convergence properties of normalized eigenfrequencies (ωa/2πc) calculated for the modes X1 (left), X2 (middle), and X3 (right) of Sample S2. The Models A, B, C, and C’ are plotted as crosses, squares, circles, and triangles, respectively.

Equations (52)

Equations on this page are rendered with MathJax. Learn more.

ɛ ( x , y ) = m , n = + ɛ mn e i ( m G x x + n G y y ) ,
y H z = i ω ɛ 0 ɛ E x ,
x H z = i ω ɛ 0 ɛ E y ,
x E y y E x = i ω μ 0 H z ,
[ y , x ] [ E x E y ] = i ω μ 0 H z ,
ɛ [ E x E y ] = 1 i ω ɛ 0 [ y x ] H z .
[ D ˜ x D ˜ y ] = ɛ [ E x E y ] = [ ɛ xx ɛ xy ɛ yx ɛ yy ] [ E x E y ] ,
( y η xx y + x η yx y + y η xy x x η yy x ) H z = ω 2 c 2 H z ,
H z ( x , y ) = m , n = + h z , mn e i ( p m x + q n y ) ,
[ D ˜ x ] = [ ɛ E x ] = [ [ ɛ ] ] [ E x ] ,
[ D ˜ y ] = [ ɛ E y ] = [ [ ɛ ] ] [ E y ] .
[ [ η ] ] A = [ [ [ ɛ ] ] 1 [ [ 0 ] ] [ [ 0 ] ] [ [ ɛ ] ] 1 ] .
[ E x E y ] = F [ E u E v ] ,
[ D ˜ u ] = [ [ 1 / ɛ ] ] 1 [ E u ] ,
[ D ˜ v ] = [ [ ɛ ] ] [ E v ] .
F = [ cos ϕ sin ϕ sin ϕ cos ϕ ] ,
[ [ E x ] [ E y ] ] = [ [ F ] ] [ [ E u ] [ E v ] ] ,
[ [ F ] ] = [ [ [ c ] ] [ [ s ] ] [ [ s ] ] [ [ c ] ] ] .
[ [ E x ] [ E y ] ] = [ [ F ] ] [ [ [ 1 ɛ ] ] [ [ 0 ] ] [ [ 0 ] ] [ [ ɛ ] ] 1 ] [ [ F 1 ] ] [ [ D ˜ x ] [ D ˜ y ] ] ,
[ [ η ] ] B = [ [ F ] ] [ [ [ 1 ɛ ] ] [ [ 0 ] ] [ [ 0 ] ] [ [ ɛ ] ] 1 ] [ [ F 1 ] ] = [ [ [ 1 ɛ ] ] [ [ c 2 ] ] + [ [ ɛ ] ] 1 [ [ s 2 ] ] , [ [ 1 ɛ ] ] [ [ c s ] ] [ [ ɛ ] ] 1 [ [ c s ] ] [ [ 1 ɛ ] ] [ [ cs ] ] [ [ ɛ ] ] 1 [ [ cs ] ] , [ [ 1 ɛ ] ] [ [ s 2 ] ] + [ [ ɛ ] ] 1 [ [ c 2 ] ] ] ,
u = [ ξ ζ ] , v = [ ζ * ξ * ]
u = e i θ [ cos θ sin θ sin θ cos θ ] [ cos i sin ] ,
θ ( r , ϕ ) = ϕ ,
( r , ϕ ) = { π 8 ( 1 + cos π r R ) ( r R ) π 8 { 1 + cos π [ r + D ( ϕ ) 2 R ] D ( ϕ ) R } ( r > R ) .
D ( ϕ ) = a / 2 max ( | cos ϕ | , | sin ϕ | )
lim r 0 u = lim r D ( ϕ ) u = 1 2 [ 1 i ] ,
[ [ η ] ] C = [ [ [ 1 ɛ ] ] [ [ ξ ξ * ] ] + [ [ ɛ ] ] 1 [ [ ζ ζ * ] ] , [ [ 1 ɛ ] ] [ [ ξ ζ * ] ] [ [ ɛ ] ] 1 [ [ ξ ζ * ] ] [ [ 1 ɛ ] ] [ [ ξ * ζ ] ] [ [ ɛ ] ] 1 [ [ ξ * ζ ] ] , [ [ 1 ɛ ] ] [ [ ζ ζ * ] ] + [ [ ɛ ] ] 1 [ [ ξ ξ * ] ] ] .
D ( ϕ ) = a / 2 max n = 0 , , 5 [ cos ( ϕ n π 3 ) ] ,
θ b ( ϕ ) = θ ( D ( ϕ ) , ϕ ) = π 2 round ( ϕ / π 2 ) ,
b ( ϕ ) = ( D ( ϕ ) , ϕ ) = π 8 ( 1 cos 4 ϕ )
θ ( r , ϕ ) = 1 2 { θ b ( ϕ ) + ϕ + [ θ b ( ϕ ) ϕ ] cos π [ r + D ( ϕ ) 2 R ] D ( ϕ ) R } ,
( r , ϕ ) = b ( ϕ ) 2 { 1 + cos π [ r + D ( ϕ ) 2 R ] D ( ϕ ) R } .
ɛ ( x , y ) = m , n = + ɛ mn e i ( m G x x + n G y y ) ,
f ( x , y ) = m , n = + f mn e i ( p m x + q n y ) .
h ( x , y ) = ɛ ( x , y ) f ( x , y ) ,
g x ( x , y ) = x f ( x , y ) ,
g y ( x , y ) = y f ( x , y ) ,
h ( x , y ) = m , n = + h mn e i ( p m x + q n y ) ,
g x ( x , y ) = m , n = + g x , mn e i ( p m x + q n y ) ,
g y ( x , y ) = m , n = + g y , mn e i ( p m x + q n y ) ,
h mn = k , l = + ɛ m k , n l f kl ,
g x , mn = i p m f mn ,
g y , mn = i q n f mn .
α ( m , n ) = m + M + 1 + ( n + N ) ( 2 M + 1 ) ,
n ( α ) = ( α 1 ) div ( 2 M + 1 ) N ,
m ( α ) = ( α 1 ) mod ( 2 M + 1 ) M ,
[ h ] = [ [ ɛ ] ] [ f ] ,
[ g x ] = i p [ f ] ,
[ g y ] = i q [ f ] ,
[ [ ɛ ] ] α β = ɛ m ( α ) m ( β ) , n ( α ) n ( β ) ,
p α β = p m ( α ) δ α β ,
q α β = q n ( α ) δ α β ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.