OSA's Digital Library

Journal of the Optical Society of America B

Journal of the Optical Society of America B

| OPTICAL PHYSICS

  • Editor: Henry Van Driel
  • Vol. 26, Iss. 11 — Nov. 1, 2009
  • pp: 1984–1993
« Show journal navigation

Dirichlet-to-Neumann map method for analyzing crossed arrays of circular cylinders

Yumao Wu and Ya Yan Lu  »View Author Affiliations


JOSA B, Vol. 26, Issue 11, pp. 1984-1993 (2009)
http://dx.doi.org/10.1364/JOSAB.26.001984


View Full Text Article

Acrobat PDF (287 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

An efficient and accurate computational method is developed for analyzing finite layers of crossed arrays of circular cylinders, including woodpile structures as special cases. The method relies on marching a few operators (approximated by matrices) from one side of the structure to another. The marching step makes use of the Dirichlet-to-Neumann (DtN) maps for two-dimensional unit cells in each layer where the structure is invariant in the direction of the cylinder axes. The DtN map is an operator that maps two wave field components to their normal derivatives on the boundary of the unit cell, and they can be easily constructed by vector cylindrical waves. Unlike existing numerical methods for crossed gratings, our method does not require a discretization of the structure. Compared with the multipole method that uses vector cylindrical wave expansions and scattering matrices, our method is relatively simple since it does not need sophisticated lattice sums techniques.

© 2009 Optical Society of America

1. INTRODUCTION

Photonic crystals (PhCs) [1

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

] have been intensively investigated both theoretically and experimentally in the last two decades. Three-dimensional (3D) PhCs with a complete bandgap in the optical wavelength region have potential applications such as ultrahigh-quality factor cavities and zero-threshold lasers. The woodpile structures [2

2. K. M. Ho, C. T. Chan, C. M. Soukoulis, R. Biswas, and M. Sigalas, “Photonic band gaps in three dimensions: new layer-by-layer periodic structures,” Solid State Commun. 89, 413–416 (1994). [CrossRef]

, 3

3. H. S. Sözüer and J. P. Dowling, “Photonic band calculations for woodpile structures,” J. Mod. Opt. 41, 231–239 (1994). [CrossRef]

, 4

4. S. Y. Lin, J. G. Fleming, D. L. Hetherington, B. K. Smith, R. Biswas, K. M. Ho, M. M. Sigalas, W. Zubrzycki, S. R. Kurtz, and J. Bur, “A three-dimensional photonic crystal operating at infrared wavelengths,” Nature 394, 251–253 (1998). [CrossRef]

] composed of alternating layers of rods have attracted much attention due to their relatively simple fabrication process compared with other 3D PhCs.

2. PROBLEM FORMULATION

We consider the diffraction of a plane wave incident upon finite layers of crossed arrays of circular cylinders. In each layer, there is a periodic array of parallel circular cylinders with period L. In a Cartesian coordinate system {x,y,z}, these cylinders are either parallel to the x axis or parallel to the y axis, depending on the layer, and they are bounded by the two planes at z=0 and z=D, where D is positive. In Fig. 1 , we show a woodpile structure as a special crossed array of cylinders where the axes of the cylinders in two layers separated by one intermediate layer are offset by half a period. Outside the cylinders and in the half-spaces z<0 and z>D, the medium is isotropic and homogeneous. We use the superscripts (1) and (2) to indicate parameters and quantities in the top (z>D) and bottom (z<0) half-spaces, respectively.

Assuming that the time dependence is exp(iωt) where ω is the angular frequency, the electromagnetic field satisfies frequency-domain Maxwell’s equations:
×E=ik0μH̃,×H̃=ik0εE,
(1)
where H̃=μ0ε0H is a scaled magnetic field, k0=ωc is the free space wavenumber, c=1ε0μ0 is the speed of light in vacuum, ɛ is the relative permittivity or dielectric constant, and μ is the relative permeability. From Eq. (1), it is straightforward to eliminate the z components and derive the following system:
z[uv]=[A11A12A21A22][uv],
(2)
where u and v are the x and y components, respectively, i.e.,
u=[ExH̃x],v=[EyH̃y],
(3)
and Apq (for p,q=1,2) are 2 × 2 matrix operators involving x and y derivatives:
A11=1ik0[0x(ε1y)x(μ1y)0],
A12=1ik0[0k02μx(ε1x)k02ε+x(μ1x)0],
A21=1ik0[0k02μ+y(ε1y)k02εy(μ1y)0],
A22=1ik0[0y(ε1x)y(μ1x)0].

Since the structure is periodic in the x and y directions and the electromagnetic field is quasi-periodic, we can formulate a boundary value problem on the square cylinder
Ω={(x,y,z)|0<x<L,0<y<L,0<z<D}.
(16)
For that purpose, we rewrite the quasi-periodic condition as
[u(L,y,z)v(L,y,z)]=eiα0L[u(0,y,z)v(0,y,z)],
(17)
[u(x,L,z)v(x,L,z)]=eiβ0L[u(x,0,z)v(x,0,z)].
(18)
Since Eq. (2) is a first-order system with respect to z, the boundary conditions at z=0 and z=D should not involve the z derivatives of u and v. To write down these boundary conditions, we define two linear operators acting on quasi-periodic functions by
B(r){ei(αjx+βky)el}=Bjk(r){ei(αjx+βky)el},
(19)
B(t){ei(αjx+βky)el}=Bjk(t){ei(αjx+βky)el},
(20)
for l=1,2 and arbitrary integers j and k, where e1 and e2 are the two columns of the 2 × 2 identity matrix, and Bjk(r) and Bjk(t) are 2 × 2 matrices given in Eqs. (14, 15). Since the operators B(r) and B(t) are linear, the superposition principle holds; therefore,
u(r)+B(r)v(r)=0,z>D,
u(t)+B(t)v(t)=0,z<0.
For z<0, we have u=u(t) and v=v(t). For z>D, we can eliminate the reflected wave and obtain
u+B(r)v=u(i)+B(r)v(i),z>D.
Since u and v are continuous across the planes at z=0 and z=D, we have the following boundary conditions:
u+B(r)v=[u(i)+B(r)v(i)]z=D+,z=D,
(21)
u+B(t)v=0,z=0.
(22)
The boundary value problem on Ω consists of Eq. (2) and boundary conditions [Eqs. (17, 18, 21, 22)].

3. OPERATOR MARCHING METHOD

Although a general numerical method such as the finite element method [14

14. G. Bao, P. Li, and H. Wu, “An adaptive edge element method with perfectly matched absorbing layers for wave scattering by biperiodic structures,” Math. Comput. 70, 1–34 (2010).

] can be used to solve the boundary value problem on Ω directly, it is possible to develop more efficient methods by utilizing the geometry features of the structure. For crossed arrays of circular cylinders, the structure can be divided into a number of layers where each layer contains a parallel array of circular cylinders. In the following, we present an operator marching (OM) method that avoids repeated calculations in different layers with identical index profiles. The main idea is to march some operators (to be approximated by matrices) from one side of the structure (z=0) to another (z=D+) and find the reflected and transmitted waves afterwards. The OM method developed in this section is an extension of the method originally developed for 2D structures [18

18. Y. Huang and Y. Y. Lu, “Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann maps,” J. Lightwave Technol. 24, 3448–3453 (2006). [CrossRef]

, 27

27. Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing interpenetrating cylinder arrays in a triangular lattice,” J. Opt. Soc. Am. B 25, 1466–1473 (2008). [CrossRef]

, 30

30. Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing periodic arrays of cylinders with oblique incident waves,” J. Opt. Soc. Am. B 26, 1442–1449 (2009). [CrossRef]

].

4. MATRIX APPROXIMATIONS FOR OPERATORS

To implement the OM method, it is necessary to approximate the operators P, Q, X, and Y by matrices. These operators act on column vectors, of which the components are quasi-periodic functions of x and y. The quasi-periodic functions are associated with α0 and β0 (the x and y components of the incident wave vector) and satisfy a condition like Eq. (7). Since the quasi-periodic functions can be expanded in Fourier series, these operators can be represented by matrices in Fourier space.

Let T be a linear operator such that if g is a vector of two quasi-periodic functions, then so is h = T g. In physical space, if one period of x or y is discretized by N points, then g and h are approximated by vectors of length 2N2, and the operator T is approximated by a (2N2)×(2N2) matrix. A different approach is to consider the action of T on different Fourier modes. For a given integer pair (j,k), we have a Fourier mode ei(αjx+βky), where αj and βk are given in Eq. (10). Since T acts on vectors, we need to consider T{ei(αjx+βky)es} for s=1,2, where e1 and e2 are the two columns of the 2 × 2 identity matrix. Since the results are also quasi-periodic, we have the following expansion:
T{ei(αjx+βky)es}=l,m=Tlmjk{ei(αlx+βmy)es}
(27)
for s=1,2, where Tlmjk are 2 × 2 matrices. Let g and h be also expanded in Fourier series as
g(x,y)=j,k=ĝjkei(αjx+βky),
h(x,y)=j,k=ĥjkei(αjx+βky),
where ĝjk and ĥjk are column vectors of length 2, then T g = h implies
j,kTlmjkĝjk=ĥlm.
(28)

If we truncate the integers j, k, l, and m to N terms (i.e., from N2 to N21 if N is even, or from (N1)2 to (N1)2 if N is odd), then all retained 2 × 2 matrices Tlmjk give a total of 4N4 numbers, and they can be stored in a (2N2)×(2N2) matrix. To actually write down this matrix, we need to order the Fourier coefficients of g and h. We introduce a column vector ĝ of length 2N2 by grouping the Fourier coefficients of g with the same k together, that is,
ĝ=[ĝ1ĝ0ĝ1],forĝk=[ĝ1,kĝ0kĝ1k].
(29)
Notice that ĝk above is a column vector of length 2N. A similar vector ĥ is defined for the Fourier coefficients of h. Based on this ordering, we can assemble the 2 × 2 matrices Tlmjk into a (2N2)×(2N2) matrix T̂ such that T g = h or Eq. (28) is approximated by T̂ĝ=ĥ. More precisely, we consider (j,k) as one index and (l,m) as another index based on the above ordering, then T̂ is a N2×N2 block matrix where each block is a 2 × 2 matrix, and Tlmjk is the block at row (l,m) and column (j,k). With such a matrix representation, operations for the operators are turned to standard matrix operations. For example, if three operators T1, T2, and T3 satisfy T1T2=T3, then the corresponding matrices satisfy T̂1T̂2=T̂3.

The differential operators Apq (p,q=1,2) in the governing Eq. (2) have simple matrix approximations in Fourier space when ɛ and μ are constants. For example, A11 acting on the Fourier mode ei(αjx+βky) gives
A11{ei(αjx+βky)es}=αjβkik0[0ε1μ10]{ei(αjx+βky)es}
(32)
for s=1,2. If we introduce 2 × 2 matrices (A11)lmjk as the general case in Eq. (27), then these matrices are nonzero only if (j,k)=(l,m), and (A11)jkjk is given in the right-hand side of Eq. (32). The matrix Â11 is assembled from the 2 × 2 matrices (A11)lmjk, and it is a block diagonal matrix with 2 × 2 blocks. The cases for Â12, Â21, and Â22 are similar; they are also block diagonal matrices with the diagonal blocks given by
(A12)jkjk=k02εμαj2ik0[0ε1μ10],
(33)
(A21)jkjk=k02εμβk2ik0[0ε1μ10],
(34)
(A22)jkjk=αjβkik0[0ε1μ10]=(A11)jkjk.
(35)
For setting up the boundary conditions [Eqs. (21, 22)], we introduced the operators B(r) and B(t) in Section 2. Their matrix representations in Fourier space are also block diagonal matrices. The diagonal blocks of B̂(r) and B̂(t) are the 2 × 2 matrices Bjk(r) and Bjk(t) given in Eqs. (14, 15), respectively.

5. THE FIRST AND LAST STEPS

In this section, we give some details for the first and last steps using the Fourier space representations of the operators. The first step is to initialize the operators Q and Y at z0 (where z0=0). From the definition, it is obvious that Y(z0) is the identity operator. Therefore, Ŷ(z0) is the identity matrix. As discussed in Section 2, for z<z0, the total wave is the transmitted wave, thus the y components of the electromagnetic field are given by
v(r)=j,k=bjk(t)exp[i(αjx+βkyγjk(2)z)]
(36)
for z<z0. Clearly, the z derivative of v is
vz(r)=ij,k=γjk(2)I2bjk(t)exp[i(αjx+βkyγjk(2)z)]
(37)
for z<z0, where I2 is the 2 × 2 identity matrix. If we follow the expansion [Eq. (27)] for operator T and define an operator S(2) based on the 2 × 2 matrices
[S(2)]lmjk={iγjk(2)I2,if(l,m)=(j,k),0,if(l,m)(j,k),}
(38)
then v satisfies
vz=S(2)v,z<z0.
(39)
The above can be applied to z=z0; therefore, we can simply set Q(z0)=S(2). In Fourier space, Q̂(z0) is a (2N2)×(2N2) matrix obtained by assembling 2 × 2 matrices [S(2)]lmjk together.

The last step is to calculate the reflected wave at zM+ (where zM=D) and the transmitted wave at z0, assuming that we know the Fourier representations of {Q, Y} or {P, X} at zM. If Q̂ and Ŷ are known at zM+, we use the y components of the electromagnetic field. For z>zM, we define an operator S(1) through its 2 × 2 matrices by replacing γjk(2) with γjk(1) in Eq. (38); then the incident and reflected waves satisfy
v(i)z=S(1)v(1),v(r)z=S(1)v(r),z>zM.
Therefore, the total field v=v(i)+v(r) satisfies
vz+S(1)v=2S(1)v(i),z>zM.
(40)
Using Q(zM+) and evaluating the above equation at z=zM+, we have
[Q(zM+)+S(1)]|v|z=zM=|2S(1)v(i)|z=zM+.
In Fourier space, the above becomes
[Q̂(zM+)+Ŝ(1)]v̂(zM)=2Ŝ(1)v̂(i)(zM+).
(41)
We can solve v̂(zM) from the above equation, then find the reflected wave by subtracting the incident wave, and find the transmitted wave from |Y(zM)v|z=zM=|v|z=z0 or
Ŷ(zM)v̂(zM)=v̂(z0).

6. TRANSITION THROUGH AN INTERFACE

In the operator marching scheme, we need to transfer the operator P or Q from zm to zm+. This is necessary if z=zm is a material interface, where ɛ or μ are discontinuous. Consider the transition process for Q. From the second equation of the first order system [Eq. (2)], we obtain
Q(zm±)v=A21(zm±)u+A22(zm±)v,
where u and v at evaluated zm, Apq(zm±) are matrix differential operators evaluated at zm±. From the equation at zm, we can solve u and insert it into the equation at zm+. This gives rise to
Q(zm+)=A21(zm+)A211(zm)[Q(zm)A22(zm)]+A22(zm+).
(42)

The above formula is especially easy to evaluate if ɛ and μ at zm± are constants, corresponding to the case where the media above and below the interfaces (in the vicinity of zm) are isotropic and homogeneous. The actions of A21(zm±) and A22(zm±) on Fourier mode ei(αjx+βky) are given in Eqs. (34, 35), if we replace ɛ and μ by ε± and μ± for their values at zm+ and zm. For the operator Gm=A21(zm+)A211(zm), then,
Gm{ei(αjx+βky)es}=k02ε+μ+βk2k02εμβk2[εε+00μμ+]{ei(αjx+βky)es}
for s=1,2. From the above, we can easily write down matrix representations of the involved operators in Fourier space, i.e., Â22(zm±) and Ĝm, then Q̂(zm+) can be evaluated by
Q̂(zm+)=Ĝm[Q̂(zm)Â22(zm)]+Â22(zm+).
(43)

7. MARCHING THROUGH A LAYER

In this section, we describe the key step where the operators are marched through a layer. Without loss of generality, we consider the first layer (z0<z<z1), which consists of a periodic array of circular cylinders parallel to the y axis. For this layer, we calculate the operators Q and Y at z1, assuming that they are given at z0. Since the electromagnetic field is quasi-periodic and the structure is y invariant, we can expand the field in Fourier series for the y variable as
v(r)=kvk(r)=kΨk(x,z)eiβky
(44)
for z0<z<z1, where r=(x,y,z). In the above, vk and a similarly defined uk satisfy the Maxwell’s equations [Eq. (2)], but they have a simple y dependence eiβky.

The OM scheme in [30

30. Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing periodic arrays of cylinders with oblique incident waves,” J. Opt. Soc. Am. B 26, 1442–1449 (2009). [CrossRef]

] was presented in the physical space based on the above Mk. In this paper, we use the Fourier space representations of the operators, because they are more convenient to carry out the switching between {P, X} and {Q, Y}. In Fourier space, Mk corresponds to a matrix M̂k satisfying
M̂k[v̂k(0)v̂k(1)]=[zv̂k(0+)zv̂k(1)],
(47)
where v̂k(0) is a column vector of the Fourier coefficients of vk(r) (as a quasi-periodic function of only x) evaluated at z0, and it is similar to the vector ĝk given in Eq. (29); the elements of zv̂k(0+) are the z derivatives of these coefficients evaluated at z0+, and v̂k(1) and zv̂k(1) are corresponding vectors at z1 or z1. When N Fourier modes are retained as before, v̂k(0) and v̂k(1) are vectors of length 2N, and M̂k is a (4N)×(4N) matrix.

From these (4N)×(4N) matrices M̂k, we obtain a (4N2)×(4N2) matrix M̂ satisfying
[zv̂(z0+)zv̂(z1)]=M̂[v̂(z0)v̂(z1)]=[M̂11M̂12M̂21M̂22][v̂(z0)v̂(z1)],
(49)
where v̂(z0) and v̂(z1) are vectors of the Fourier coefficients for v at z0 and z1, respectively, and the blocks M̂11,M̂12, are themselves block diagonal matrices given by
M̂pq=[M̂1,pqM̂0,pqM̂1,pq]
(50)
for p,q=1,2.

The definitions of Q and Y given in Eqs. (25, 26) imply that
Q̂(zm±)v̂(zm)=zv̂(zm±),Ŷ(zm)v̂(zm)=v̂(z0).
Therefore, Eq. (49) can be replaced by
[Q̂(z0+)00Q̂(z1)][v̂(z0)v̂(z1)]=[M̂11M̂12M̂21M̂22][v̂(z0)v̂(z1)].
Solving v̂(z0) from the first equation above, we can easily derive the following formulas:
Q̂(z1)=M̂22+M̂21[Q̂(z0+)M̂11]1M̂12,
(51)
Ŷ(z1)=Ŷ(z0)[Q̂(z0+)M̂11]1M̂12.
(52)

For woodpile structures, the axes of the cylinders in two layers separated by one intermediate layer are offset by a distance of L2 in the direction perpendicular to the cylinder axes. Therefore, the 2D unit cell
Ω3={(x,z)|0<x<L,z2<z<z3}
for the third layer contains two half-cylinders with their centers located at x=0 and x=L, respectively. For such a layer, we consider shifted unit cell
Ω3={(x,z)|L2<x<3L2,z2<z<z3}.
From the DtN map Λk for Ω3, we calculate the reduced DtN map Mk by eliminating the vertical edges. After that, we can use a shifting technique to find the reduced DtN map Mk for Ω3 [26

26. Y. Huang and Y. Y. Lu, “Modeling photonic crystals with complex unit cells by Dirichlet-to-Neumann maps,” J. Comput. Math. 25, 337–349 (2007).

, 27

27. Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing interpenetrating cylinder arrays in a triangular lattice,” J. Opt. Soc. Am. B 25, 1466–1473 (2008). [CrossRef]

, 30

30. Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing periodic arrays of cylinders with oblique incident waves,” J. Opt. Soc. Am. B 26, 1442–1449 (2009). [CrossRef]

]. Alternatively, we can find M̂k directly from Mk using a discrete Fourier transform based on the shifted points:
xl=L2+xl=L2+(l0.5)LN,1lN.

If the layer consists of an array of cylinders parallel to the x axis, we need to use the operators P, X and the vector u. In that case, the second ordering of the Fourier coefficients given in Eq. (30) is preferred; therefore, the marching step is carried out with P̃, X̃ and ũ.

8. SWITCHING THE OPERATORS

Between layers of different orientations, it is necessary to switch between the operator pairs {Q, Y} and {P, X}. From the definitions of P and Q and the governing Eq. (2), we have
[P00Q][uv]=[A11A12A21A22][uv].
Similar to the derivation of Eq. (51), we first solve v in terms of u as
v=(QA22)1A21u,
(53)
and then obtain
P=A11+A12(QA22)1A21.
(54)
In Fourier space and following the first ordering [Eq. (29)], the above becomes
P̂=Â11+Â12(Q̂Â22)1Â21.
(55)
The matrices P̂,Q̂,Â11,Â12, are all (2N2)×(2N2) matrices where N is the number of retained terms for Fourier series in x or y. The switching formula [Eq.(55)] is only used at the interfaces between the layers where the medium is homogeneous, i.e., ɛ and μ are constants. In that case, Â11,Â12, are quite simple, as discussed in Section 4.

Since the marching step for a x invariant layer uses P̃ and W̃ (or X̃) based on the second ordering [Eq. (30)] of the Fourier coefficients, it is necessary to transform P̂ to P̃, etc. This can easily be done by a permutation process given in Eq. (31). After marching through a x invariant layer, it is necessary to transform Q̃ and W̃ back to P̂ and Ŷ. The process is similar to the one given above.

9. NUMERICAL EXAMPLES

To illustrate our method, we consider some examples based on the parameters used by Yasumoto and Jia in [17

17. K. Yasumoto and H. Jia, “Electromagnetic scattering from multilayered crossed-arrays of circular cylinders,” SPIE5445, 200–205 (2004).

]. The structure consists of M layers of crossed arrays of circular dielectric cylinders surrounded by air. The dielectric constant and the radius of the cylinders are ɛ = 5 and r=0.25L, respectively, where L is the period in both x and y directions. The distance between two nearby layers is also L, that is, zmzm1=L for 1mM. Since the structure is nonmagnetic, we have μ = 1 everywhere. As in Section 3, we assume that the cylinders in adjacent layers are orthogonal to each other. More precisely, the cylinders in the even and odd layers are parallel to the x and y axes, respectively. We consider two cases. The first case is a regular crossed array of cylinders where all even (or odd) layers are identical. Therefore, the structure repeats itself every two layers. The second case is a woodpile structure, where the axes of cylinders in two layers separated by one intermediate layer are offset by L2. Therefore, the woodpile structure repeats itself every four layers. For both cases we calculate the reflection and transmission spectra for plane waves at normal incidence given in the top region z>D=zM. Therefore, the wave vector of the incident wave is (α0,β0,γ00(1))=(0,0,k0). The incident wave is normalized and given by
a(i)=[10],b(i)=[01],
where a(i) and b(i) are the coefficients of the incident wave as given in Eq. (4). For ωL(2πc)1 and the normal incident wave, only the (0,0)th diffraction order is propagating. Therefore, we are mainly interested in the coefficients a00(r), b00(r), a00(t), and b00(t), as defined in Eqs. (8, 9). In the following, we show the power reflectance and transmittance of the (0,0)th diffraction order given by
R00=a00(r)2=b00(r)2,T00=a00(t)2=b00(t)2,
where ‖⋅‖ denotes the magnitude of a vector.

In Fig. 2 , we show the transmission spectra for M=4. The top and bottom plots are results for a regular crossed arrays structure and a woodpile structure, respectively. Although there are only four layers, the transmission spectra already exhibit rapid variations at some frequencies. In Fig. 3 , we show the reflection spectra for M=32. Once again, the top and bottom plots are the reflection spectra for a regular crossed array structure and a woodpile structure, respectively. We can easily identify a number of frequency intervals where the reflectance is nearly 100%. Outside these intervals, the reflection spectra exhibit rapid oscillations as frequency varies and they tend to be more oscillatory for larger frequencies. The results for the two structures (regular crossed arrays and woodpile) are certainly different, and the difference is more pronounced for larger frequencies. But overall, the difference is not as significant as we have expected.

The case for the regular crossed arrays and M=32 was previously analyzed by Yasumoto and Jia [17

17. K. Yasumoto and H. Jia, “Electromagnetic scattering from multilayered crossed-arrays of circular cylinders,” SPIE5445, 200–205 (2004).

] using a multipole method (also called cylindrical harmonic expansion method) that involves scattering matrix and lattice sums techniques. The results shown in the top plot of Fig. 3 roughly agree with those reported in [17

17. K. Yasumoto and H. Jia, “Electromagnetic scattering from multilayered crossed-arrays of circular cylinders,” SPIE5445, 200–205 (2004).

]. In particular, there are good agreement for the high-reflectance intervals. However, the rapid oscillatory behavior of the spectra was not revealed in [17

17. K. Yasumoto and H. Jia, “Electromagnetic scattering from multilayered crossed-arrays of circular cylinders,” SPIE5445, 200–205 (2004).

]. This may be caused by a coarse sampling of the frequency.

To validate our results, we check the power balance law: R00+T00=1 for ωL(2πc)1. The numerical results in Figs. 2, 3 are obtained with N=9 or N=10. For these results, the difference between R00+T00 and 1 is on the order of 107 or 105 for low or high frequencies, respectively. We also check the numerical convergence for increasing N at fixed frequencies. Since exact solutions are not available, we use the results obtained with a large N as a reference solution for computing approximate relative errors. For R00, the approximate relative error is defined as EN=|R00(N)R00(20)|R00(20) for N15. In Fig. 4 , we show the approximate relative errors for four different cases. It is observed that the errors appear to decrease exponentially as N is increased.

10. CONCLUSIONS

In this paper, we developed a Dirichlet-to-Neumann (DtN) operator marching (OM) method for analyzing crossed arrays of circular cylinders including woodpile structures as special cases. For a 3D multilayer structure where each layer is a periodic array of circular cylinders, we used cylindrical waves to construct the DtN maps of the 2D unit cells for each Fourier mode (from an expansion in the variable along the axes of the cylinders) and then applied these DtN maps in an OM scheme to calculate transmission and reflection spectra. The OM scheme involves some operators that are approximated by (2N2)×(2N2) matrices in Fourier space, where N can be quite small. Typically, three significant digits can be obtained using N=9. Since analytic solutions are used to construct the DtN maps, a discretization of the structure is not needed. Furthermore, the same DtN map applies to all layers with the same index profile. Therefore, our method is efficient for multilayer structures that have some periodicity in the z direction. For simplicity, the method is presented for the special case where the periods in the x and y directions are the same, but it can be easily extended to the case where the periods in these two directions are different. Compared with the multipole method developed in [16

16. G. H. Smith, L. C. Botten, R. C. McPhedran, and N. A. Nicorovici, “Cylinder gratings in conical incidence with applications to woodpile structures,” Phys. Rev. E 67, 056620 (2003). [CrossRef]

, 17

17. K. Yasumoto and H. Jia, “Electromagnetic scattering from multilayered crossed-arrays of circular cylinders,” SPIE5445, 200–205 (2004).

], our method is relatively simple, since we do not require sophisticated lattice sums techniques.

ACKNOWLEDGMENTS

This research was partially supported by a grant from the Research Grants Council of Hong Kong Special Administrative Region, China, project CityU 102207.

Fig. 1 Woodpile structure as special crossed arrays of cylinders.
Fig. 2 Transmission spectra of a four-layer regular crossed array structure (top plot) and a four-layer woodpile structure (bottom plot).
Fig. 3 Reflection spectra of a 32-layer regular crossed array structure (top plot) and a 32-layer woodpile structure (bottom plot).
Fig. 4 Relative errors of the power reflectance R00 for regular crossed arrays (two left plots) and woodpile structures (two right plots).
1.

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

2.

K. M. Ho, C. T. Chan, C. M. Soukoulis, R. Biswas, and M. Sigalas, “Photonic band gaps in three dimensions: new layer-by-layer periodic structures,” Solid State Commun. 89, 413–416 (1994). [CrossRef]

3.

H. S. Sözüer and J. P. Dowling, “Photonic band calculations for woodpile structures,” J. Mod. Opt. 41, 231–239 (1994). [CrossRef]

4.

S. Y. Lin, J. G. Fleming, D. L. Hetherington, B. K. Smith, R. Biswas, K. M. Ho, M. M. Sigalas, W. Zubrzycki, S. R. Kurtz, and J. Bur, “A three-dimensional photonic crystal operating at infrared wavelengths,” Nature 394, 251–253 (1998). [CrossRef]

5.

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]

6.

S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8, 173–190 (2001). [CrossRef] [PubMed]

7.

D. C. Dobson, J. Gopalakrishnan, and J. E. Pasciak, “An efficient method for band structure calculations in 3D photonic crystals,” J. Comput. Phys. 161, 668–679 (2000). [CrossRef]

8.

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-difference Time Domain Method, 2nd ed. (Artech House, 2000).

9.

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

10.

E. Popov and M. Nevière, “Maxwell equations in Fourier space: a fast-converging formulation for diffraction by arbitrary shaped, periodic, anisotropic media,” J. Opt. Soc. Am. A 18, 2886–2894 (2001). [CrossRef]

11.

M. Nevière and E. Popov, Light Propagation in Periodic Media (Marcel Dekker, 2003).

12.

J. Elschner, R. Hinder, and G. Schmidt, “Finite element solution of conical diffraction problems,” Adv. Comput. Math. 16, 139–156 (2002). [CrossRef]

13.

G. Bao, Z. M. Chen, and H. J. Wu, “Adaptive finite-element method for diffraction gratings,” J. Opt. Soc. Am. A 22, 1106–1114 (2005). [CrossRef]

14.

G. Bao, P. Li, and H. Wu, “An adaptive edge element method with perfectly matched absorbing layers for wave scattering by biperiodic structures,” Math. Comput. 70, 1–34 (2010).

15.

E. Popov, M. Nevière, B. Gralak, and G. Tayeb, “Staircase approximation validity for arbitrary-shaped gratings,” J. Opt. Soc. Am. A 19, 33–42 (2002). [CrossRef]

16.

G. H. Smith, L. C. Botten, R. C. McPhedran, and N. A. Nicorovici, “Cylinder gratings in conical incidence with applications to woodpile structures,” Phys. Rev. E 67, 056620 (2003). [CrossRef]

17.

K. Yasumoto and H. Jia, “Electromagnetic scattering from multilayered crossed-arrays of circular cylinders,” SPIE5445, 200–205 (2004).

18.

Y. Huang and Y. Y. Lu, “Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann maps,” J. Lightwave Technol. 24, 3448–3453 (2006). [CrossRef]

19.

J. Yuan and Y. Y. Lu, “Photonic bandgap calculations using Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. A 23, 3217–3222 (2006). [CrossRef]

20.

S. Li and Y. Y. Lu, “Multipole Dirichlet-to-Neumann map method for photonic crystals with complex unit cells,” J. Opt. Soc. Am. A 24, 2438–2442 (2007). [CrossRef]

21.

J. Yuan, Y. Y. Lu, and X. Antoine, “Modeling photonic crystals by boundary integral equation and Dirichlet-to-Neumann maps,” J. Comput. Phys. 9, 4617–4629 (2008). [CrossRef]

22.

H. Xie and Y. Y. Lu, “Modeling two-dimensional anisotropic photonic crystals by Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. A 26, 1606–1614 (2009). [CrossRef]

23.

J. Yuan and Y. Y. Lu, “Computing photonic band structures by Dirichlet-to-Neumann maps: the triangular lattice,” Opt. Commun. 273, 114–120 (2007). [CrossRef]

24.

Y. Huang, Y. Y. Lu, and S. Li, “Analyzing photonic crystal waveguides by Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. B 24, 2860–2867 (2007). [CrossRef]

25.

S. Li and Y. Y. Lu, “Computing photonic crystal defect modes by Dirichlet-to-Neumann maps,” Opt. Express 15, 14454–14466 (2007). [CrossRef] [PubMed]

26.

Y. Huang and Y. Y. Lu, “Modeling photonic crystals with complex unit cells by Dirichlet-to-Neumann maps,” J. Comput. Math. 25, 337–349 (2007).

27.

Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing interpenetrating cylinder arrays in a triangular lattice,” J. Opt. Soc. Am. B 25, 1466–1473 (2008). [CrossRef]

28.

Z. Hu and Y. Y. Lu, “Efficient analysis of photonic crystal devices by Dirichlet-to-Neumann maps,” Opt. Express 16, 17383–17399 (2008). [CrossRef] [PubMed]

29.

Z. Hu and Y. Y. Lu, “Improved Dirichlet-to-Neumann map method for modeling extended photonic crystal devices,” Opt. Quantum Electron. 40, 921–932 (2008). [CrossRef]

30.

Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing periodic arrays of cylinders with oblique incident waves,” J. Opt. Soc. Am. B 26, 1442–1449 (2009). [CrossRef]

OCIS Codes
(000.4430) General : Numerical approximation and analysis
(050.5298) Diffraction and gratings : Photonic crystals

ToC Category:
Diffraction and Gratings

History
Original Manuscript: July 10, 2009
Manuscript Accepted: August 24, 2009
Published: October 5, 2009

Virtual Issues
October 8, 2009 Spotlight on Optics

Citation
Yumao Wu and Ya Yan Lu, "Dirichlet-to-Neumann map method for analyzing crossed arrays of circular cylinders," J. Opt. Soc. Am. B 26, 1984-1993 (2009)
http://www.opticsinfobase.org/josab/abstract.cfm?URI=josab-26-11-1984


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals: Modeling the Flow of Light (Princeton Univ. Press, 1995).
  2. K. M. Ho, C. T. Chan, C. M. Soukoulis, R. Biswas, and M. Sigalas, “Photonic band gaps in three dimensions: new layer-by-layer periodic structures,” Solid State Commun. 89, 413-416 (1994). [CrossRef]
  3. H. S. Sözüer and J. P. Dowling, “Photonic band calculations for woodpile structures,” J. Mod. Opt. 41, 231-239 (1994). [CrossRef]
  4. S. Y. Lin, J. G. Fleming, D. L. Hetherington, B. K. Smith, R. Biswas, K. M. Ho, M. M. Sigalas, W. Zubrzycki, S. R. Kurtz, and J. Bur, “A three-dimensional photonic crystal operating at infrared wavelengths,” Nature 394, 251-253 (1998). [CrossRef]
  5. 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]
  6. S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell's equations in a planewave basis,” Opt. Express 8, 173-190 (2001). [CrossRef] [PubMed]
  7. D. C. Dobson, J. Gopalakrishnan, and J. E. Pasciak, “An efficient method for band structure calculations in 3D photonic crystals,” J. Comput. Phys. 161, 668-679 (2000). [CrossRef]
  8. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-difference Time Domain Method, 2nd ed. (Artech House, 2000).
  9. L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc. Am. A 14, 2758-2767 (1997). [CrossRef]
  10. E. Popov and M. Nevière, “Maxwell equations in Fourier space: a fast-converging formulation for diffraction by arbitrary shaped, periodic, anisotropic media,” J. Opt. Soc. Am. A 18, 2886-2894 (2001). [CrossRef]
  11. M. Nevière and E. Popov, Light Propagation in Periodic Media (Marcel Dekker, 2003).
  12. J. Elschner, R. Hinder, and G. Schmidt, “Finite element solution of conical diffraction problems,” Adv. Comput. Math. 16, 139-156 (2002). [CrossRef]
  13. G. Bao, Z. M. Chen, and H. J. Wu, “Adaptive finite-element method for diffraction gratings,” J. Opt. Soc. Am. A 22, 1106-1114 (2005). [CrossRef]
  14. G. Bao, P. Li, and H. Wu, “An adaptive edge element method with perfectly matched absorbing layers for wave scattering by biperiodic structures,” Math. Comput. 70, 1-34 (2010).
  15. E. Popov, M. Nevière, B. Gralak, and G. Tayeb, “Staircase approximation validity for arbitrary-shaped gratings,” J. Opt. Soc. Am. A 19, 33-42 (2002). [CrossRef]
  16. G. H. Smith, L. C. Botten, R. C. McPhedran, and N. A. Nicorovici, “Cylinder gratings in conical incidence with applications to woodpile structures,” Phys. Rev. E 67, 056620 (2003). [CrossRef]
  17. K. Yasumoto and H. Jia, “Electromagnetic scattering from multilayered crossed-arrays of circular cylinders,” SPIE 5445, 200-205 (2004).
  18. Y. Huang and Y. Y. Lu, “Scattering from periodic arrays of cylinders by Dirichlet-to-Neumann maps,” J. Lightwave Technol. 24, 3448-3453 (2006). [CrossRef]
  19. J. Yuan and Y. Y. Lu, “Photonic bandgap calculations using Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. A 23, 3217-3222 (2006). [CrossRef]
  20. S. Li and Y. Y. Lu, “Multipole Dirichlet-to-Neumann map method for photonic crystals with complex unit cells,” J. Opt. Soc. Am. A 24, 2438-2442 (2007). [CrossRef]
  21. J. Yuan, Y. Y. Lu, and X. Antoine, “Modeling photonic crystals by boundary integral equation and Dirichlet-to-Neumann maps,” J. Comput. Phys. 9, 4617-4629 (2008). [CrossRef]
  22. H. Xie and Y. Y. Lu, “Modeling two-dimensional anisotropic photonic crystals by Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. A 26, 1606-1614 (2009). [CrossRef]
  23. J. Yuan and Y. Y. Lu, “Computing photonic band structures by Dirichlet-to-Neumann maps: the triangular lattice,” Opt. Commun. 273, 114-120 (2007). [CrossRef]
  24. Y. Huang, Y. Y. Lu and S. Li, “Analyzing photonic crystal waveguides by Dirichlet-to-Neumann maps,” J. Opt. Soc. Am. B 24, 2860-2867 (2007). [CrossRef]
  25. S. Li and Y. Y. Lu, “Computing photonic crystal defect modes by Dirichlet-to-Neumann maps,” Opt. Express 15, 14454-14466 (2007). [CrossRef] [PubMed]
  26. Y. Huang and Y. Y. Lu, “Modeling photonic crystals with complex unit cells by Dirichlet-to-Neumann maps,” J. Comput. Math. 25, 337-349 (2007).
  27. Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing interpenetrating cylinder arrays in a triangular lattice,” J. Opt. Soc. Am. B 25, 1466-1473 (2008). [CrossRef]
  28. Z. Hu and Y. Y. Lu, “Efficient analysis of photonic crystal devices by Dirichlet-to-Neumann maps,” Opt. Express 16, 17383-17399 (2008). [CrossRef] [PubMed]
  29. Z. Hu and Y. Y. Lu, “Improved Dirichlet-to-Neumann map method for modeling extended photonic crystal devices,” Opt. Quantum Electron. 40, 921-932 (2008). [CrossRef]
  30. Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing periodic arrays of cylinders with oblique incident waves,” J. Opt. Soc. Am. B 26, 1442-1449 (2009). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1 Fig. 2 Fig. 3
 
Fig. 4
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited