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

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

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,

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

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

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:

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 $2N^2$, and the operator T is approximated by a $\left(2N^2\right)\times \left(2N^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:

If we truncate the integers j, k, l, and m to N terms (i.e., from $-N∕2$ to $N∕2-1$ if N is even, or from $-(N-1)∕2$ to $(N-1)∕2$ if N is odd), then all retained 2 × 2 matrices $T_{lm}^{jk}$ give a total of $4N^4$ numbers, and they can be stored in a $\left(2N^2\right)\times \left(2N^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 $2N^2$ by grouping the Fourier coefficients of g with the same k together, that is,

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

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

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

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

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,

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

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

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

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