## 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] 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, 3, 4] composed of alternating layers of rods have attracted much attention due to their relatively simple fabrication process compared with other 3D PhCs.

To analyze PhCs and related devices, numerical methods are indispensable. Band structures of perfectly periodic PhCs, as well as waveguides and microcavities in PhCs with defects, give rise to eigenvalue problems that can be solved by many existing methods such as the plane wave expansion method [5, 6] and the finite element method [7]. PhC devices such as waveguide bends, branches, and couplers involving waveguide and cavities give rise to more challenging boundary-value problems. In particular, for PhCs that are finite in one direction, it is important to calculate the transmission and reflection spectra for given incident waves. The finite-difference time domain (FDTD) method [8] is widely used to study these problems, but its accuracy and efficiency are often limited. Notice that the FDTD requires a small mesh size to resolve curved material interfaces (especially when the index contrast is high), and it often has difficulties truncating periodic structures that extend to infinity. The mathematical problem associated with a finite PhC, which is still infinite and periodic in the other directions, is identical to the problem for diffraction gratings. A woodpile structure consisting of a finite number of layers gives rise to the same 3D mathematical problem as crossed gratings. Therefore, existing numerical methods for diffraction gratings, such as the Fourier modal method (FMM) [9], the differential method [10, 11], and the finite element method [12, 13, 14] can be used to analyze finite woodpile structures. However, these methods are typically quite expensive to use. In particular, if the woodpile is composed of circular rods, FMM and other modal methods must approximate the structure by many uniform layers, but the involved staircase approximation may have convergence problems [15]. The differential method gives rise to a boundary value problem for a large system of ordinary differential equations, and it is still expensive to solve. While the finite element method is extremely general, it gives rise to large, complex, non-Hermitian, indefinite, and sparse linear systems that are difficult to solve. The multipole method [16, 17] is more efficient because it uses vector cylindrical waves to take advantage of the special geometric features. However, since each layer has infinite number of cylinders, the multipole method requires sophisticated lattice sums techniques.

Recently, an efficient computational method for PhCs and related devices was developed based on the so-called Dirichlet-to-Neumann (DtN) maps of unit cells [18, 19]. The DtN map of a unit cell is an operator that maps the wave field to its normal derivative on the boundary of the unit cell, and it can be approximated by a small matrix based on special solutions obtained analytically or numerically [20, 21, 22]. Using the DtN maps, we can significantly reduce the number of unknowns by avoiding the interiors of unit cells. For ideal two-dimensional (2D) PhCs composed of parallel and infinitely long cylinders and for waves propagating in the plane perpendicular to the cylinder axes, the DtN map method has been used to calculate band structures [19, 23], waveguide modes [24], cavity modes [25], transmission and reflection spectra of finite PhCs [18, 26, 27], and to analyze general PhC devices connected by a few PhC waveguides [28, 29]. In a previous work [30], we extended the DtN map method to out-of-plane propagation problems of 2D PhCs. Using vector cylindrical waves, we constructed matrix DtN maps for two field components and calculated transmission and reflection spectra for oblique plane incident waves. In this paper, we extend the DtN map method to 3D structures with a finite number of layers of crossed arrays of circular cylinders, including woodpile structures as special cases.

## 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 $\mathrm{exp}(-i\omega t)$ where ω is the angular frequency, the electromagnetic field satisfies frequency-domain Maxwell’s equations:

*z*components and derive the following system:

**u**and

**v**are the

*x*and

*y*components, respectively, i.e.,

*x*and

*y*derivatives:

For $z>D$, we specify a plane incident wave as

*x*and

*y*can be observed from the inverse of ${B}_{00}^{\left(i\right)}$:

In connection with the plane incident wave and the periodicity of the structure in the *x* and *y* directions, the electromagnetic field is quasi-periodic in *x* and *y*; that is,

**Φ**and

**Ψ**are periodic in

*x*and

*y*with period

*L*. For $z>D$, the total electromagnetic field is the sum of the incident and reflected waves. For $z<0$, the total field is the transmitted wave. It is well known that the reflected and transmitted waves can be expanded as

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

*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

*j*and

*k*, where ${\mathbf{e}}_{1}$ and ${\mathbf{e}}_{2}$ are the two columns of the 2 × 2 identity matrix, and ${B}_{jk}^{\left(r\right)}$ and ${B}_{jk}^{\left(t\right)}$ are 2 × 2 matrices given in Eqs. (14, 15). Since the operators ${\mathcal{B}}^{\left(r\right)}$ and ${\mathcal{B}}^{\left(t\right)}$ are linear, the superposition principle holds; therefore,

**u**and

**v**are continuous across the planes at $z=0$ and $z=D$, we have the following boundary conditions:

## 3. OPERATOR MARCHING METHOD

Although a general numerical method such as the finite element method [14] 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, 27, 30].

Consider a multilayer structure where the layers are separated by $0={z}_{0}<{z}_{1}<\dots <{z}_{M}=D$, and the $m\text{th}$ layer is specified by ${z}_{m-1}<z<{z}_{m}$. Each layer is invariant in one direction (*x* or *y*) and periodic in the other direction (*y* or *x*) with a period *L*. At ${z}_{m}$ or ${z}_{m}^{\pm}$ (for $0\le m\le M$), we define four operators *P*, *Q*, *X* and *Y,*:satisfying

**u**and

**v**satisfying Eq. (2), the quasi-periodic conditions [Eqs. (17, 18)], and the boundary condition [Eq. (22)]. If $z={z}_{m}$ is a material interface,

**u**and

**v**are continuous at ${z}_{m}$, but their

*z*derivatives are not, so therefore

*X*and

*Y*can be defined at ${z}_{m}$ but

*P*and

*Q*are defined at both ${z}_{m}^{-}$ and ${z}_{m}^{+}$. To understand the above definitions, we notice that Eq. (2) is a first-order system with respect to

*z*. Therefore, in general, if

**u**is specified at ${z}_{m}$, it should have a unique solution satisfying Eqs. (17, 18, 22). For that solution, we can evaluate ${\partial}_{z}\mathbf{u}$ at ${z}_{m}^{\pm}$ and

**u**at ${z}_{0}$, and they are related to

**u**at ${z}_{m}$ by the linear operators $\mathcal{P}\left({z}_{m}^{\pm}\right)$ and $\mathcal{X}\left({z}_{m}\right)$. The four operators are not independent. In fact, the operators

*Q*and

*Y*are related to

*P*and

*X*, respectively. However, it is convenient to work with the

*x*(or

*y*) components of the electromagnetic field in layers where the cylinders are parallel to the

*x*(or

*y*) axis. Therefore, we introduce these four operators and use them alternatively in layers with different axial orientations.

In the following, we assume that the cylinders in the $m\text{th}$ layer are parallel to the *x* axis if *m* is an even integer and parallel to the *y* axis if *m* is odd. Then, our OM method proceeds with the following steps:

*M*. The four operators above are approximated by $\left(2{N}^{2}\right)\times \left(2{N}^{2}\right)$ matrices in Fourier space, where

*N*is the number of retained modes in Fourier series for functions of

*x*or

*y*. The details are given in Section 4. The first and last steps are related to the definitions [Eqs. (23, 24, 25, 26)] and the boundary conditions at ${z}_{0}^{-}$ and ${z}_{M}^{+}$, and they are explained in Section 5. Steps 2 and 3 involve three different tasks: transition through an interface, marching through a layer, and switching the operators. We explain these three tasks in Sections 6, 7, 8, respectively. To march a pair of operators through a layer, we use DtN maps for unit cells of this layer. This is the key step that allows us to take full advantage of the geometric features.

## 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 ${\alpha}_{0}$ and ${\beta}_{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 $2{N}^{2}$, and the operator *T* is approximated by a $\left(2{N}^{2}\right)\times \left(2{N}^{2}\right)$ 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 ${e}^{i({\alpha}_{j}x+{\beta}_{k}y)}$, where ${\alpha}_{j}$ and ${\beta}_{k}$ are given in Eq. (10). Since *T* acts on vectors, we need to consider $\mathcal{T}\left\{{e}^{i({\alpha}_{j}x+{\beta}_{k}y)}{\mathbf{e}}_{s}\right\}$ for $s=1,2$, where ${\mathbf{e}}_{1}$ and ${\mathbf{e}}_{2}$ are the two columns of the 2 × 2 identity matrix. Since the results are also quasi-periodic, we have the following expansion:

**g**and

**h**be also expanded in Fourier series as

*T*

**g**=

**h**implies

If we truncate the integers *j*, *k*, *l*, and *m* to *N* terms (i.e., from $-N\u22152$ to $N\u22152-1$ if *N* is even, or from $-(N-1)\u22152$ to $(N-1)\u22152$ if *N* is odd), then all retained 2 × 2 matrices ${T}_{lm}^{jk}$ give a total of $4{N}^{4}$ numbers, and they can be stored in a $\left(2{N}^{2}\right)\times \left(2{N}^{2}\right)$ matrix. To actually write down this matrix, we need to order the Fourier coefficients of **g** and **h**. We introduce a column vector $\widehat{\mathbf{g}}$ of length $2{N}^{2}$ by grouping the Fourier coefficients of **g** with the same *k* together, that is,

**h**. Based on this ordering, we can assemble the 2 × 2 matrices ${T}_{lm}^{jk}$ into a $\left(2{N}^{2}\right)\times \left(2{N}^{2}\right)$ matrix $\widehat{\mathcal{T}}$ such that

*T*

**g**=

**h**or Eq. (28) is approximated by $\widehat{\mathcal{T}}\phantom{\rule{0.1em}{0ex}}\widehat{\mathbf{g}}=\widehat{\mathbf{h}}$. More precisely, we consider $(j,k)$ as one index and $(l,m)$ as another index based on the above ordering, then $\widehat{\mathcal{T}}$ is a ${N}^{2}\times {N}^{2}$ block matrix where each block is a 2 × 2 matrix, and ${T}_{lm}^{jk}$ 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 ${\mathcal{T}}_{1}$, ${\mathcal{T}}_{2}$, and ${\mathcal{T}}_{3}$ satisfy ${\mathcal{T}}_{1}{\mathcal{T}}_{2}={\mathcal{T}}_{3}$, then the corresponding matrices satisfy ${\widehat{\mathcal{T}}}_{1}{\widehat{\mathcal{T}}}_{2}={\widehat{\mathcal{T}}}_{3}$.

We can also order the Fourier coefficients by grouping those with the same *j* together. For the coefficients of **g**, we let

*S*satisfying $S={S}^{-1}$. We have

The differential operators ${\mathcal{A}}_{pq}$ $(p,q=1,2)$ in the governing Eq. (2) have simple matrix approximations in Fourier space when ɛ and μ are constants. For example, ${\mathcal{A}}_{11}$ acting on the Fourier mode ${e}^{i({\alpha}_{j}x+{\beta}_{k}y)}$ gives

## 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 ${z}_{0}^{-}$ (where ${z}_{0}=0$). From the definition, it is obvious that $\mathcal{Y}\left({z}_{0}^{-}\right)$ is the identity operator. Therefore, $\widehat{\mathcal{Y}}\left({z}_{0}^{-}\right)$ is the identity matrix. As discussed in Section 2, for $z<{z}_{0}$, the total wave is the transmitted wave, thus the *y* components of the electromagnetic field are given by

*z*derivative of

**v**is

*T*and define an operator ${\mathcal{S}}^{\left(2\right)}$ based on the 2 × 2 matrices

**v**satisfies

The last step is to calculate the reflected wave at ${z}_{M}^{+}$ (where ${z}_{M}=D$) and the transmitted wave at ${z}_{0}^{-}$, assuming that we know the Fourier representations of {*Q*, *Y*} or {*P*, *X*} at ${z}_{M}$. If $\widehat{\mathcal{Q}}$ and $\widehat{\mathcal{Y}}$ are known at ${z}_{M}^{+}$, we use the *y* components of the electromagnetic field. For $z>{z}_{M}$, we define an operator ${\mathcal{S}}^{\left(1\right)}$ through its 2 × 2 matrices by replacing ${\gamma}_{jk}^{\left(2\right)}$ with ${\gamma}_{jk}^{\left(1\right)}$ in Eq. (38); then the incident and reflected waves satisfy

## 6. TRANSITION THROUGH AN INTERFACE

In the operator marching scheme, we need to transfer the operator *P* or *Q* from ${z}_{m}^{-}$ to ${z}_{m}^{+}$. This is necessary if $z={z}_{m}$ 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

**u**and

**v**at evaluated ${z}_{m}$, ${\mathcal{A}}_{pq}\left({z}_{m}^{\pm}\right)$ are matrix differential operators evaluated at ${z}_{m}^{\pm}$. From the equation at ${z}_{m}^{-}$, we can solve

**u**and insert it into the equation at ${z}_{m}^{+}$. This gives rise to

The above formula is especially easy to evaluate if ɛ and μ at ${z}_{m}^{\pm}$ are constants, corresponding to the case where the media above and below the interfaces (in the vicinity of ${z}_{m}$) are isotropic and homogeneous. The actions of ${\mathcal{A}}_{21}\left({z}_{m}^{\pm}\right)$ and ${\mathcal{A}}_{22}\left({z}_{m}^{\pm}\right)$ on Fourier mode ${e}^{i({\alpha}_{j}x+{\beta}_{k}y)}$ are given in Eqs. (34, 35), if we replace ɛ and *μ* by ${\epsilon}^{\pm}$ and ${\mu}^{\pm}$ for their values at ${z}_{m}^{+}$ and ${z}_{m}^{-}$. For the operator ${\mathcal{G}}_{m}={\mathcal{A}}_{21}\left({z}_{m}^{+}\right){\mathcal{A}}_{21}^{-1}\left({z}_{m}^{-}\right)$, then,

## 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 $({z}_{0}<z<{z}_{1})$, 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 ${z}_{1}$, assuming that they are given at ${z}_{0}$. 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

*y*dependence ${e}^{i{\beta}_{k}y}$.

For such a *y*-invariant layer and a fixed integer *k*, the DtN map method developed in our previous work [30] is applicable. Consider the 2D unit cell of the first layer,

*x*direction, i.e., ${\mathbf{v}}_{k}(x+L,y,z)={e}^{i{\alpha}_{0}L}{\mathbf{v}}_{k}(x,y,z)$, we can eliminate ${\mathbf{v}}_{k}$ and its

*x*derivative on the vertical edges and find the reduced DtN map ${\mathcal{M}}_{k}$, satisfying

*x*by

*N*points as ${x}_{l}=(l-0.5)L\u2215N$ for $1\le l\le N$, then ${\mathbf{v}}_{k}^{\left(0\right)}$ denotes a column vector of length $2N$ consisting of ${\mathbf{v}}_{k}({x}_{l},y,{z}_{0})$ for $1\le l\le N$, ${\partial}_{z}{\mathbf{v}}_{k}^{(0+)}$ denotes a vector for the corresponding

*z*derivatives at ${z}_{0}^{+}$, and ${\mathbf{v}}_{k}^{\left(1\right)}$ and ${\partial}_{z}{\mathbf{v}}_{k}^{(1-)}$ are similarly defined. Notice that ${\mathcal{M}}_{k}$ is represented by a $\left(4N\right)\times \left(4N\right)$ matrix. More details can be found in [30].

The OM scheme in [30] was presented in the physical space based on the above ${\mathcal{M}}_{k}$. 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, ${\mathcal{M}}_{k}$ corresponds to a matrix ${\widehat{\mathcal{M}}}_{k}$ satisfying

*x*) evaluated at ${z}_{0}$, and it is similar to the vector ${\widehat{\mathbf{g}}}_{k}$ given in Eq. (29); the elements of ${\partial}_{z}{\widehat{\mathbf{v}}}_{k}^{(0+)}$ are the

*z*derivatives of these coefficients evaluated at ${z}_{0}^{+}$, and ${\widehat{\mathbf{v}}}_{k}^{\left(1\right)}$ and ${\partial}_{z}{\widehat{\mathbf{v}}}_{k}^{(1-)}$ are corresponding vectors at ${z}_{1}$ or ${z}_{1}^{-}$. When

*N*Fourier modes are retained as before, ${\widehat{\mathbf{v}}}_{k}^{\left(0\right)}$ and ${\widehat{\mathbf{v}}}_{k}^{\left(1\right)}$ are vectors of length $2N$, and ${\widehat{\mathcal{M}}}_{k}$ is a $\left(4N\right)\times \left(4N\right)$ matrix.

From ${\mathcal{M}}_{k}$ we can calculate ${\widehat{\mathcal{M}}}_{k}$ by a simple extension of the discrete Fourier transform. If *g* is a scalar quasi-periodic function of *x* and it is evaluated at ${x}_{l}$ for $1\le l\le N$, then we can approximate its Fourier coefficients by the discrete Fourier transform

*N*terms as before. The two column vectors for $g\left({x}_{l}\right)$ and ${\widehat{g}}_{j}$, respectively, are linked by an $N\times N$ matrix with entries ${e}^{i{\alpha}_{j}{x}_{l}}$. Now for vector quasi-periodic function ${\mathbf{v}}_{k}$, its values at the

*N*discrete points of

*x*(and at fixed

*y*and

*z*) are related to the Fourier coefficients given in the vector ${\widehat{\mathbf{v}}}_{k}$ by a $\left(2N\right)\times \left(2N\right)$ matrix

*F*, which is also an $N\times N$ block matrix with 2 × 2 blocks. The $(l,j)$ block of

*F*is ${F}_{lj}={e}^{i{\alpha}_{j}{x}_{l}}{I}_{2}$, where ${I}_{2}$ is the 2 × 2 identity matrix. Let ${\mathcal{M}}_{k}$ and ${\widehat{\mathcal{M}}}_{k}$ be given in 2 × 2 blocks as

From these $\left(4N\right)\times \left(4N\right)$ matrices ${\widehat{\mathcal{M}}}_{k}$, we obtain a $\left(4{N}^{2}\right)\times \left(4{N}^{2}\right)$ matrix $\widehat{\mathcal{M}}$ satisfying

**v**at ${z}_{0}$ and ${z}_{1}$, respectively, and the blocks ${\widehat{\mathcal{M}}}_{11},{\widehat{\mathcal{M}}}_{12},\dots $ are themselves block diagonal matrices given by

The definitions of *Q* and *Y* given in Eqs. (25, 26) imply that

For woodpile structures, the axes of the cylinders in two layers separated by one intermediate layer are offset by a distance of $L\u22152$ in the direction perpendicular to the cylinder axes. Therefore, the 2D unit cell

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 cellFrom the DtN map ${\Lambda}_{k}^{\prime}$ for ${\Omega}_{3}^{\prime}$, we calculate the reduced DtN map ${\mathcal{M}}_{k}^{\prime}$ by eliminating the vertical edges. After that, we can use a shifting technique to find the reduced DtN map ${\mathcal{M}}_{k}$ for ${\Omega}_{3}$ [26, 27, 30]. Alternatively, we can find ${\widehat{\mathcal{M}}}_{k}$ directly from ${\mathcal{M}}_{k}^{\prime}$ using a discrete Fourier transform based on the shifted points: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 $\stackrel{\u0303}{\mathcal{P}}$, $\stackrel{\u0303}{\mathcal{X}}$ and $\stackrel{\mathrm{\u0303}}{\mathbf{u}}$.

## 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

**v**in terms of

**u**asand then obtain

*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, ${\widehat{\mathcal{A}}}_{11},{\widehat{\mathcal{A}}}_{12},\dots $ are quite simple, as discussed in Section 4.

The operators *X* and *Y* are only needed for computing the transmitted waves. While we can switch between *Y* and *X*, it is advantageous to use a different operator *W*, so that the switching between $\text{\hspace{0.17em}}{\phantom{|}\mathbf{u}|}_{z={z}_{0}}$ and $\text{\hspace{0.17em}}{\phantom{|}\mathbf{v}|}_{z={z}_{0}}$ can be avoided. We define *W* at ${z}_{m}$ by

**v**is related to

**u**by Eq. (53), we haveIn Fourier space, the above formula becomes

Since the marching step for a *x* invariant layer uses $\stackrel{\u0303}{\mathcal{P}}$ and $\stackrel{\u0303}{\mathcal{W}}$ (or $\stackrel{\u0303}{\mathcal{X}}$) based on the second ordering [Eq. (30)] of the Fourier coefficients, it is necessary to transform $\widehat{\mathcal{P}}$ to $\stackrel{\u0303}{\mathcal{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 $\stackrel{\u0303}{\mathcal{Q}}$ and $\stackrel{\u0303}{\mathcal{W}}$ back to $\widehat{\mathcal{P}}$ and $\widehat{\mathcal{Y}}$. 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]. 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, ${z}_{m}-{z}_{m-1}=L$ for $1\le m\le M$. 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 $L\u22152$. 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={z}_{M}$. Therefore, the wave vector of the incident wave is $({\alpha}_{0},{\beta}_{0},-{\gamma}_{00}^{\left(1\right)})=(0,0,-{k}_{0})$. The incident wave is normalized and given by

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] 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]. In particular, there are good agreement for the high-reflectance intervals. However, the rapid oscillatory behavior of the spectra was not revealed in [17]. This may be caused by a coarse sampling of the frequency.

To validate our results, we check the power balance law: ${R}_{00}+{T}_{00}=1$ for $\omega L\u2215\left(2\pi c\right)\le 1$. The numerical results in Figs. 2, 3 are obtained with $N=9$ or $N=10$. For these results, the difference between ${R}_{00}+{T}_{00}$ and 1 is on the order of ${10}^{-7}$ or ${10}^{-5}$ 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 ${R}_{00}$, the approximate relative error is defined as ${\mathrm{E}}_{N}=|{R}_{00}^{\left(N\right)}-{R}_{00}^{\left(20\right)}|\u2215{R}_{00}^{\left(20\right)}$ for $N\le 15$. 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 $\left(2{N}^{2}\right)\times \left(2{N}^{2}\right)$ 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, 17], 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.

**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]