## Abstract

An efficient numerical method based on the Dirichlet-to-Neumann (DtN) maps of the unit cells is developed for accurate simulations of two-dimensional photonic crystal (PhC) devices in the frequency domain. The DtN map of a unit cell is an operator that maps the wave field on the boundary of the cell to its normal derivative and it can be approximated by a small matrix. Using the DtN maps of the regular and defect unit cells, we can avoid computations in the interiors of the unit cells and calculate the wave field only on the edges. This gives rise to a significant reduction in the total number of unknowns. Reasonably accurate solutions can be obtained using 10 to 15 unknowns for each unit cell. In contrast, standard finite element, finite difference or plane wave expansion methods may require a few hundreds unknowns for each unit cell at the same level of accuracy. We illustrate our method by a number of examples, including waveguide bends, branches, microcavities coupled with waveguides, waveguides with stubs, etc.

©2008 Optical Society of America

## 1. Introduction

In recent years, photonic crystals (PhCs) [1–3] have been extensively studied both theoretically and experimentally, due to their unusual ability to control and manipulate light. Because of the periodicity of the dielectric constant, PhCs exhibit unusual dispersion properties and frequency gaps in which propagating Bloch waves do not exist. These properties have been widely used to design photonic crystal devices, such as waveguide bends [4–6], branches [7], frequency filters [8], waveguide couplers [7], Mach-Zehnder interferometers [9], etc. Numerical methods are essential to analyze basic properties of PhCs and to design and optimize PhC devices. Fundamental problems such as band structures, waveguide and cavity modes lead to eigenvalue problems that can be solved by a variety of different numerical methods [10–24]. PhC devices, such as a PhC waveguide bend, usually give rise to more challenging boundary value problems. While the band structure problem is formulated on a unit cell, a PhC device has to be studied on a much larger domain with more complicated boundary conditions.

Many PhC devices are simulated in time domain, for example, by the finite difference time domain (FDTD) method. For some problems, such as the propagation of a pulse, the time domain methods are essential. Other problems, such as the transmission and reflection spectra, are more naturally formulated in the frequency domain. However, even for two-dimensional (2D) problems, standard numerical methods for frequency domain formulations, such as the finite element method, often give rise to large linear systems that are complex, non-Hermitian, indefinite but sparse. These systems are expensive to solve by direct methods. Iterative methods often have a very slow convergence and even fail to converge, since existing preconditioning techniques for indefinite systems are not very effective. However, frequency domain formulations have one important advantage: they allow us to take advantage of the geometric features of the structure.

Consider a 2D PhC composed of a lattice of infinitely long and parallel cylinders in a homogeneous background, such as air-holes in a dielectric medium or dielectric rods in air. When cavities and waveguides are introduced as point and line defects, the structure loses its periodicity, but it still has many identical unit cells. Very often, there are only two different types of unit cells: the regular unit cell and the defect unit cell. Of course, the wave fields are different on different cells, but it is possible to take advantage of the many identical cells by using their Dirichlet-to-Neumann (DtN) maps.

For a given domain Ω and a linear homogeneous differential equation for some function *u*, the DtN map Λ is an operator that maps *u* on the boundary of Ω to the normal derivative of *u* on the same boundary. For PhCs, the domain Ω is chosen as a unit cell. With the DtN maps of the unit cells, we can write down equations for the wave field on the edges of the unit cells only. Therefore, the DtN maps allow us to avoid solving the wave field in the interiors of the unit cells completely. In previous works, the DtN maps of the unit cells have been used to develop efficient methods for computing band structures [21,22], waveguide modes [23], cavity modes [24] and transmission/reflection spectra [25–27] of finite PhCs. In this paper, the DtN-map method is extended to general boundary value problems for arbitrary 2D PhC devices in an infinite background PhC. The device is allowed to have a finite number of PhC waveguides that extend to infinity. Interfaces between PhC and non-PhC structures are not considered in this paper. We restrict our attention to pure 2D structures that are invariant in the third direction and assume that the waves are propagating in a 2D plane. The problems associated with PhC slabs are certainly very important, but they are not studied here.

## 2. Equations on cell edges

For pure 2D structures which are invariant in the *z* direction and for waves propagating in the *xy* plane, the governing equation is the Helmholtz equation

where *k*
_{0} is the free space wavenumber, *n*=*n*(x) is the refractive index function and *x*=(*x*,*y*). For the *E*-polarization, *u* is the *z*-component of the electric field and *ρ*=1. For the *H*-polarization, *u* is the *z*-component of the magnetic field and *ρ*=*n*
^{2}. We consider 2D PhC structures that can be divided into many unit cells. Usually, there are only a small number of distinct unit cells corresponding to the original bulk PhC and defects. Let Ω be a unit cell and Γ be its boundary, we first find the Dirichlet-to-Neumann (DtN) map Λ so that

where *u* satisfies the Helmholtz Eq. (1) in Ω and *ν* is a unit normal vector of Γ. The operator Λ can be approximated by a *K*×*K* matrix, if we choose *K* sampling points on Γ. As described in [25] and [21], we can find the matrix approximation of Λ by assuming that the general solution of Eq. (1) in Ω can be approximated by a sum of *K* special solutions. That is

where *ϕ*
* _{j}* satisfies Eq. (1) in Ω. If we evaluate

*ϕ*

*and the normal derivative of*

_{j}*φ*

*at the*

_{j}*K*points on Γ, we can eliminate the unknown coefficients {

*c*

*} and find the matrix approximation of Λ. For unit cells containing circular cylinders, we choose*

_{j}*ϕ*

*as a cylindrical wave which can be written down analytically. For more complicated unit cells, the DtN map can be approximated by the methods developed in [28, 29].*

_{j}With the DtN map Λ, we can write down an equation for each edge of the unit cells. To do this, we need to be more specific in choosing the points and the unit normal vector on Γ. For a square unit cell Ω_{1} given by 0<*x*,*y*<*a*, if we choose *N* points on each edge, we can order *u* at the 4*N* points on Γ as a column vector. In MATLAB notation, we have

*u*=[*u*
_{01}; *v*
_{01}; *u*
_{11}; *v*
_{11}],

where *u*
_{01}, *u*
_{11}, *v*
_{01} and *v*
_{11} are column vectors of length *N* representing *u* evaluated on the four edges. More specifically, we have

${u}_{01}=\left[\begin{array}{c}u({\tau}_{1},0)\\ u({\tau}_{2},0)\\ \vdots \\ u({\tau}_{N},0)\end{array}\right],{u}_{11}=\left[\begin{array}{c}u({\tau}_{1},a)\\ u({\tau}_{2},a)\\ \vdots \\ u({\tau}_{N},a)\end{array}\right],{v}_{01}=\left[\begin{array}{c}u(0,{\tau}_{1})\\ u(0,{\tau}_{2})\\ \vdots \\ u(0,{\tau}_{N})\end{array}\right],{v}_{11}=\left[\begin{array}{c}u(a,{\tau}_{1})\\ u(a,{\tau}_{2})\\ \vdots \\ u(a,{\tau}_{N})\end{array}\right],$

where *τ*
* _{j}*=(

*j*-0.5)

*a*/

*N*for 1≤

*j*≤

*N*. For square unit cells, we also choose the unit normal vector such that

*∂*

_{ν}

*u*becomes

*∂*

_{x}

*u*and

*∂*

_{y}

*u*on the vertical and horizontal edges, respectively. Therefore, the DtN map Λ

^{(1)}gives

In the above, we have partitioned Λ^{(1)} as 4×4 blocks where each block is an *N*×*N* matrix. Similarly, for the unit cell Ω_{2} given by *a*<*x*<2*a* and 0<*y*<*a*, we have a DtN map Λ^{(2)} satisfying

${\mathrm{\Lambda}}^{\left(2\right)}\left[\begin{array}{c}{u}_{02}\\ {v}_{11}\\ {u}_{12}\\ {v}_{21}\end{array}\right]=\left[\begin{array}{c}{\partial}_{y}{u}_{02}\\ {\partial}_{x}{v}_{11}\\ {\partial}_{y}{u}_{12}\\ {\partial}_{x}{v}_{21}\end{array}\right]$,

where *u*
_{02}, *u*
_{12} and *v*
_{21} are column vectors of length *N* for *u* evaluated on the edges of Ω_{2}. On the common edge of Ω_{1} and Ω_{2} at *x*=*a*, we can evaluate *∂*
_{x}
*u*, denoted by *∂*
_{x}
*v*
_{11} in the discrete form, by Λ^{(1)} and Λ^{(2)} in Ω_{1} and Ω_{2}, respectively. The continuity of *ρ*
^{-1}
*∂*
_{x}
*u* gives rise to

$$=\frac{1}{{\rho}_{2}}\left({\mathrm{\Lambda}}_{21}^{\left(2\right)}{u}_{02}+{\mathrm{\Lambda}}_{22}^{\left(2\right)}{v}_{11}+{\mathrm{\Lambda}}_{23}^{\left(2\right)}{u}_{12}+{\mathrm{\Lambda}}_{24}^{\left(2\right)}{v}_{21}\right),$$

where *ρ* is defined in connection with Eq. (1), *ρ*
_{1} and *ρ*
_{2} are *ρ* evaluated at the left and right sides of the common edge at *x*=*a* and Λ^{(2)}
* _{jk}* are the blocks of Λ

^{(2)}. For the

*E*polarization, we have

*ρ*

_{1}=

*ρ*

_{2}=1. For the

*H*polarization, if the background media in the two unit cells are identical, we still have

*ρ*

_{1}=

*ρ*

_{2}. Eq. (5), actually a system of

*N*equations, is identified as the equation for the edge associated with

*v*

_{11}and it links the seven edges of the two neighboring unit cells. Clearly, for any interior edge, which is a common edge of two neighboring unit cells in the computation domain, we can establish a similar equation using the DtN maps of the unit cells.

## 3. Boundary conditions

We consider PhC devices developed in an infinite 2D bulk PhC. Away from a finite domain, we have a few PhC waveguides that extend to infinity. These waveguides serve as the ports where light can propagate to or away from the device. For practical numerical simulations, we have to truncate the domain and use appropriate boundary conditions. We consider boundary value problems at a given frequency which is inside a band gap of the background PhC.

If the bulk PhC is composed of a square lattice of cylinders (dielectric rods or air-holes) in a homogeneous background medium, we truncate the domain to a rectangle or a union of a few rectangles, following the edges of unit cells (both regular and defect unit cells). The boundary of the computation domain is thus composed of segments of straight lines. If such a line segment does not cut through a waveguide and it is sufficiently far away from other defect structures of the device, we can simply use a zero Dirichlet boundary condition, since we assumed that the frequency is within a band gap of the PhC. In the following, we describe boundary conditions on a line segment that cuts though one or more PhC waveguides. To simplify the presentation, let us shift the *x* axis so that the half plane *y*>0 represents a PhC waveguide whose axis is parallel to the *y* axis. It is possible that the waveguide is in fact a super-waveguide composed of a few line defects parallel to the *y* axis. The non-trivial part of the PhC device is in the lower half plane given by *y*<0. A segment of the *x*-axis is on the boundary of the computation domain. Let us also shift the *y* axis, so that the segment is given by 0<*x*<ma, where *a* is the lattice constant of the bulk PhC and *m* is an integer. The waveguide in the upper half plane is periodic in the *y* direction. For a line defect, the period is the lattice constant *a*. In Fig. 2, we show one period (in the *y* direction) of a simple line defect waveguide where the *x* variable is truncated to cover *m*=7 lattice constants. In the following, we establish a boundary condition at *y*=0 for 0<*x*<*ma* based on the Bloch modes of the waveguide in the upper half plane.

As in [23], we calculate the Bloch modes using the DtN map of the supercell which covers one period of the waveguide, i.e., 0<*y*<*a*. The DtN map *M* of the supercell satisfies

where *u*
_{0} and *u*
_{1} represent *u* evaluated at *y*=0 and *y*=*a*, respectively. Consistent with the domain truncation, the PhC waveguide in the upper half plane is truncated to 0<*x*<*ma*, assuming that the field is zero at *x*=0 and *x*=*ma*. Therefore, the supercell is composed of *m* unit cells given by

${\mathrm{\Omega}}_{j}=\{(x,y)\mid \left(j-1\right)a<x<\mathit{ja},0<y<a\}$

for 1≤*j*≤*m*. In the discrete case, with *N* points on each edge of the unit cells, *u*
_{0} and *u*
_{1} are column vectors of length *mN* corresponding to *x*
* _{k}*=(

*k*-0.5)

*a*/

*N*for 1≤

*k*≤

*mN*, and

*M*is a (2

*mN*)×(2

*mN*) matrix. Using the DtN map Λ

^{(j)}of the unit cell Ω

_{j}for 1≤

*j*≤

*m*, we can calculate the matrix

*M*by eliminating the field on the vertical edges, i.e.,

*v*

_{j}_{1}for 1≤

*j*≤

*m*. Since we assume that the field is zero at

*x*=0 and

*x*=

*ma*, we have

*v*

_{01}=

*v*

_{m1}=0. If we write down the equations for

*m*-1 interior vertical edges, such as Eq. (5) for

*v*

_{11}, we obtain the system

where *A*
_{1} is an (*m*-1)*N*×(*m*-1)*N* block tridiagonal matrix, *A*
_{2} is an (*m*-1)*N*×(2*mN*) matrix and

${v}_{1}=\left[\begin{array}{c}{v}_{11}\\ {v}_{21}\\ \vdots \\ {v}_{m-1,1}\end{array}\right],{u}_{0}=\left[\begin{array}{c}{u}_{01}\\ {u}_{02}\\ {u}_{03}\\ \vdots \\ {u}_{0m}\end{array}\right],{u}_{1}=\left[\begin{array}{c}{u}_{11}\\ {u}_{12}\\ {u}_{13}\\ \vdots \\ {u}_{1m}\end{array}\right].$.

On the other hand, for each unit cell, as it is clear from Eq. (4) for Ω_{1}, we can evaluate *∂*
_{y}
*u* on the horizontal edges using the 1st and 3rd block rows of its DtN map. This leads to

where *B*
_{1} is a (2*mN*)×(2*mN*) square matrix and *B*
_{2} is a (2*mN*)×(*m*-1)*N* matrix. If we solve *v*
_{1} from (7) and insert it into (8), we obtain Eq. (6) and the DtN map of the supercell:

The PhC waveguide in the upper half plane is periodic in *y* with period *a*. The wave field in the waveguide is a superposition of Bloch modes. A Bloch mode is a special solution of the Helmholtz equation given by

where Φ is periodic in *y* with period *a*. It is known that the Bloch modes appear in pairs. Corresponding to the Bloch mode above with the propagation constant *β*(which may be complex), there is another Bloch mode with the propagation constant -*β*, namely, *w*̃(*x*,*y*)=$\tilde{\Phi}$
(*x*,*y*)*e*
^{-iβy}, where $\tilde{\Phi}$
is periodic in *y* with period *a*. This is true without any symmetry assumptions on the waveguide and it is valid for lossy media where the refractive index may be complex. If the waveguide (periodically extended to -∞<*y*<∞) has a reflection symmetry in *y*, i.e., *n*(*x*,-*y*)=*n*(*x*,*y*), then we have $\tilde{\Phi}$
(*x*,*y*)=Φ(*x*,*-y*). For lossless media, the complex conjugates of *w* and *w*̃ are also solutions of Eq. (1), therefore, ±*β* are also propagation constants of Bloch modes. Notice that for *µ*=*e*
* ^{iβ a}*, we have

${w}_{1}=\mu {w}_{0},{\partial}_{y}{w}_{1}=\mu {\partial}_{y}{w}_{0}$,

where *w*
_{0}=*w*(*x*,0), *∂*
_{y}
*w*
_{0}=*∂*
_{y}
*w*(*x*,0), etc. Since Eq. (6) is valid for any solution of the Helmholtz equation, if we write *M* in a 2×2 block form, we can re-write Eq. (6) for the Bloch mode solution w as the following eigenvalue problem:

where *I* is the identity operator,*M*
* _{jk}* (for

*j*,

*k*=1, 2) are the blocks of

*M*and

*µ*is the eigenvalue. In the discrete case,

*w*

_{0}and

*∂*

_{y}

*w*

_{0}are column vectors of length

*mN*,

*M*

*and*

_{jk}*I*are (

*mN*)×(

*mN*) matrices.

Since the Bloch modes appear in pairs, we can choose *β*
* _{j}* such that its imaginary part is positive if it is complex, and its average power flux on the

*x*axis is positive (towards

*y*=+∞) if

*β*

*is real. Therefore, the wave field in the waveguide can be decomposed as*

_{j}*u*=

*u*

^{+}+

*u*

^{-}, where

In the above, *u*
^{+} is the outgoing (towards *y*=+∞) wave field component including propagating Bloch modes (real *β*
* _{j}*) having a positive net power flux and evanescent Bloch modes (complex

*β*

*) that decay exponentially as*

_{j}*y*is increased,

*u*

^{-}is the opposite incoming wave field component. Usually,

*u*

^{-}is assumed to be given and it may contain a single propagating Bloch mode, and

*u*

^{+}is the unknown outgoing wave field. If we evaluate

*u*

^{+}at

*y*=0 and

*y*=

*a*, we have

${u}_{0}^{+}=\sum _{j=1}^{\infty}{c}_{j}{\varphi}_{j},\phantom{\rule{.2em}{0ex}}{u}_{1}^{+}=\sum _{j=1}^{\infty}{c}_{j}{\mu}_{j}{\varphi}_{j},$

where *ϕ*
* _{j}*=Φ

*(*

_{j}*x*,0) and

*µ*

*=*

_{j}*e*

*. Let us define a linear operator*

^{iβja}*T*satisfying

$T{\varphi}_{j}={\mu}_{j}{\varphi}_{j},\phantom{\rule{.2em}{0ex}}j=1,2,\dots $

then *u*
^{+}
_{1}=*Tu*
^{+}
_{0} from the linearity of *T*. From the DtN map *M* of the supercell given in Eq. (6), we have *∂*
_{y}
*u*
^{+}
_{0}=*M*
_{11}
*u*
^{+}
_{0}+*M*
_{12}
*u*
^{+}
_{1}. This gives rise to a boundary condition for *u*
^{+}. That is

Similarly, we can derive a boundary condition for *u*
^{-} at *y*=0. For that purpose, we assume that the waveguide in the upper half plane is periodically extended to the lower half plane, and apply the DtN map *M* to *u*
^{-} on the supercell given by -*a*<*y*<0. With a linear operator *T*̃ defined by

*T*̃*ϕ*̃* _{j}*=

*µ*

_{j}*ϕ*̃

*,*

_{j}*j*=1,2, …

where *ϕ*̃* _{j}*=

*Φ*̃

*(*

_{j}*x*,0), we obtain

If the waveguide has the reflection symmetry in *y*, then *ϕ*̃* _{j}*=

*ϕ*

*, thus*

_{j}*T*̃=

*T*. From the decomposition

*u*=

*u*

^{+}+

*u*

^{-}and the conditions (13) and (14), we can eliminate

*u*

^{+}and obtain the following boundary condition for the total field:

The boundary condition is given at *y*=0^{+} to allow possible material discontinuities at *y*=0. Material interfaces could exist between the periodic waveguide in the upper half plane and the structure in the lower half plane.

In the discrete case, the eigenvalue problem (11) has 2*mN* eigenvalues µ* _{j}* and 1/µ

*for*

_{j}*j*=1, 2, …,

*mN*. The first half of the eigenvector, i.e.

*w*

_{0}in (11), corresponds to

*ϕ*

*(*

_{j}*x*) and

*ϕ*̃

*(*

_{j}*x*) for

*x*evaluated at

*x*

*for 1≤*

_{k}*k*≤

*mN*. Under the assumption that these vectors are linearly independent, we have the following explicit formula for the matrix

*T*,

$T=\left[{\varphi}_{1}{\varphi}_{2}\dots \right]\left[\begin{array}{ccc}{\mu}_{1}& \phantom{\rule{.2em}{0ex}}& \phantom{\rule{.2em}{0ex}}\\ \phantom{\rule{.2em}{0ex}}& {\mu}_{2}& \phantom{\rule{.2em}{0ex}}\\ \phantom{\rule{.2em}{0ex}}& \phantom{\rule{.2em}{0ex}}& \ddots \end{array}\right]{\left[{\varphi}_{1}{\varphi}_{2}\dots \right]}^{-1},$

where *ϕ*
* _{j}* denotes the column vector for

*ϕ*

*(*

_{j}*x*) at the discrete points. The case for

*T*̃ is similar. Finally, 𝓛

^{+}and 𝓛

^{-}are (

*mN*)×(

*mN*) matrices given as in (13) and (14).

Similarly, if the PhC waveguide is given in the lower half plane (*y*<0), the boundary condition is

where *u*
^{+} is the given incident wave in the waveguide. The boundary conditions (15) and (16) allow us to set up equations for edges on those boundary segments that terminate semi-infinite PhC waveguides. On such an edge, we can evaluate the normal derivative of the field from the interior and exterior of the computation domain separately. Within the computation domain, the edge belongs to a unit cell and the normal derivative can be evaluated using the DtN map of that unit cell. Outside the computation domain, the normal derivative is evaluated by a boundary condition such as (15) or (16), or similar ones in the x direction. An equation for this edge is then established by the continuity of *ρ*
^{-1}
*∂*
_{ν}
*u*. If the edge belongs to a boundary segment of length *ma*, where *a* is the lattice constant of the background PhC, then this equation establishes a link between the total of *m*+3 edges on the boundary segment and the unit cell. Overall, we can establish and solve a linear system for the wave field on all interior edges of the computation domain and on all boundary segments that terminate semi-infinite waveguides. The coefficient matrix is somewhat sparse, since the equation for an interior edge involves only seven edges of two neighboring unit cells.

## 4. Triangular lattice

In this section, we consider structures on 2D PhCs composed of circular cylinders in a triangular lattice. For greater symmetry, we use hexagon unit cells. The DtN map of such a cell can be constructed based on the same approximation (3) with *K*=6*N*, where *K* is the number of special solutions and *N* is the number of sampling points on each edge. To actually write down a matrix approximation for the DtN map Λ, we need to order the six edges, order the sampling points on each edge and choose a unit normal vector *ν* for each edge. For example, we can order the edges in the clockwise direction and order the points on each edge so that their *x*-coordinates are always increasing (and their *y*-coordinates are increasing if their *x*-coordinates are constant). Notice that the sampling points on opposite edges are ordered in the same direction. This allows us to avoid vector reversing when we set up equations on the edges. Finally, we also choose the same unit vector for opposite edges of the hexagon. In fact, we choose the unit normal vector as an upward vector with a positive *y* component (or as the unit vector in the positive *x* direction on a vertical edge). The DtN map of the unit cell is then written in 6×6 blocks, where each block is an *N*×*N* matrix.

Similar to the case for square lattice PhC devices, we consider only ideal structures designed on an infinite background PhC. Away from a finite domain, we have a few PhC waveguides that extend to infinity. A finite computation domain is obtained by truncating the original unbounded domain following the edges of certain hexagon unit cells. On each edge in the interior of the computation domain, we can evaluate the normal derivative of the wave field based on the DtN maps of its two neighboring unit cells. An equation for this edge is then established from the continuity of *ρ*
^{-1}
*∂*
_{ν}*u* as in section 3. This equations connects the 11 edges of two neighboring unit cells.

The boundary of the computation domain can be divided into a few curves. Some of these curves are sufficiently far away from the defect structures, so that a zero Dirichlet boundary condition can be used there. Other curves are used to terminate semi-infinite PhC waveguides that extend to infinity. A boundary condition is then needed on such a curve. Consider a semi-infinite PhC waveguide along the positive y direction. One period (in the y direction) of a single line defect waveguide is shown in Fig. 3. The transverse direction of the waveguide is also truncated. In Fig. 3, five layers of cylinders are retained on each side of the line defect. The supercell of the truncated PhC is bounded by the curves Γ_{0} (lower red curve), Γ_{1} (upper red curve), two lateral edges in the left and two lateral edges in the right. Notice that Γ_{1} is a vertical translation of Γ_{0} by the period of the waveguide *a* (also the lattice constant of the background PhC).

We can derive a boundary condition on Γ_{0} following the same steps as in section 3. First, we calculate the DtN map *M* of the supercell. If we let ${u}_{0}=u{\mid}_{{\mathrm{\Gamma}}_{0}}$
${u}_{1}=u{\mid}_{{\mathrm{\Gamma}}_{1}}$, ${\partial}_{\nu}{u}_{0}={\partial}_{\nu}u{\mid}_{{\mathrm{\Gamma}}_{0}}$ and ${\partial}_{\nu}{u}_{1}={\partial}_{\nu}u{\mid}_{{\mathrm{\Gamma}}_{1}}$, then Eq. (6) for *M* is still valid if we replace the y derivative by the normal derivative. As before, *M* is obtained from the DtN maps of the unit cells by eliminating the interior edges inside the supercell and imposing the zero boundary condition on the four lateral edges. Next, we calculate the Bloch modes of the waveguide. The eigenvalue problem is given as (11) for *w*
_{0}=*w*|_{Γ0}, but *∂*
_{y}
*w*
_{0} should be replaced by *∂*
_{ν}*w*|_{Γ0}. As before, the eigenvalues appear in pairs: µ* _{j}* and 1/µ

*for*

_{j}*j*=1, 2, 3, …, and the first half of the corresponding eigenvectors are

*ϕ*

*and*

_{j}*ϕ*̃

*, respectively. Finally, we define the operators*

_{j}*T*,

*T*̃, 𝓛

^{+}and 𝓛

^{-}as before, then the boundary condition for terminating the semi-infinite PhC waveguide in the positive

*y*direction is

where *u*
^{-} is the given incident field in the waveguide (coming down from *y*=+∞). Similarly, the boundary condition for terminating a PhC waveguide in the negative *y* direction is

where *u*
^{+} is the given incident field coming from *y*=-∞.

As in section 3, for each edge on a boundary curve that terminates a semi-infinite waveguide, we can establish an equation using one of the boundary conditions above and the DtN map of a related unit cell in the computation domain. Together with the equations for interior edges, we have a complete linear system for the wave field on edges of the unit cells in the computation domain.

## 5. Numerical examples

In this section, we illustrate our method by some numerical examples. We start with the waveguide bend proposed by Mekis *et al* [4]. The background PhC is a square lattice of dielectric rods in air. The refractive index and the radius of the rods are 3.4 and 0.18*a*, respectively, where *a* is the lattice constant. For the *E* polarization, the bulk PhC has a band gap given by 0.302<ω*a*/(2*πc*)<0.443. A straight PhC waveguide is formed by removing one row of rods. The waveguide supports one propagating mode for 0.312<ω*a*/(2*πc*)<0.443. The bend proposed in [4] is shown in Fig. 4 (left panel). The objective is to calculate the transmission and reflection properties of the bend for an incoming propagating mode in the horizontalwaveguide. Both the waveguides and the bulk PhC are assumed to extend to infinity. In the FDTD simulations of Mekis *et al* [4], a large computation domain covering 100×120=12000 unit cells was used. In our calculations, it is only necessary to use 11×11=121 unit cells (precisely as shown in Fig. 4), since the boundary conditions at the left and top edges accurately simulate the PhC waveguides that extend to infinity. On the bottom and right edges, we use a zero Dirichlet boundary condition, since the wave field decays exponentially away from the defects if the frequency is in the band gap. For this structure, there are only two distinct unit cells: the regular cell with a rod inside and the empty defect cell. The DtN maps of these two unit cells are approximated by (4*N*)×(4*N*) matrices, where *N* is the number of sampling points on each edge. This is based on approximating the general solution inside each unit cell as a sum of 4*N* cylindrical waves. Using these two DtN maps, we can construct the boundary conditions at the left and top boundaries. The condition at the left boundary is like (16), where *y* should be replaced by *x* and *u*
^{+} represents the given incoming propagating mode in the waveguide. The condition at the top boundary is like (15) with *u*
^{-}=0, since there are only outgoing waves in the vertical waveguide. In these boundary conditions, the operators L^{±} are approximated by (11*N*)×(11*N*) matrices. With the DtN maps of the unit cells and the boundary conditions, we can then set up a linear system of equations for the wave field on all edges of the unit cells. Since *N* points are used on on each edge, the total number of unknowns is 242*N*. This corresponds to 2*N* unknowns for each unit cell. Although a square unit cell has four edges, each interior edge is shared by two unit cells. In our calculations, accurate results are obtained with *N*=5. This implies that the wave field in each unit cell is represented by 4*N*=20 cylindrical waves. The results are satisfactory, since the size of the unit cell is less than one half of the free space wavelength. Furthermore, the linear system of 1210 unknowns has a sparse coefficient matrix, since the equation for each edge involves only 6 additional edges of the two neighboring unit cells. In Fig. 4 (right panel), we show the transmission and reflection spectra for the frequency range 0.32≤ω*a*/(2*πc*)≤0.44. The solid lines and the small circles in Fig. 4 are solutions obtained with *N*=5 and *N*=7, respectively, and they are indistinguishable from each other. In fact, our numerical results indicate an exponential convergence as *N* is increased. The original FDTD results in [4] have some ripples in the low frequency region. More accurate solutions were obtained by Koshiba *et al* [7] using a finite element time domain method and by Smajic *et al* [5] using a multiple multipole method. Our results are in excellent agreement with those reported in [7] and [5]. When the wave field on the boundary of a unit cell is known, we can easily calculate the field everywhere in the cell by its cylindrical wave expansion (of 4*N* terms). In Fig. 5, we show the electric field patterns for ω*a*/(2*πc*)=0.353 where near 100% transmission is observed, and for ω*a*/(2*πc*)=0.42 where the transmission is relatively low. Overall, our method is very efficient, since the number of unknowns is quite small. For comparison, the finite element time domain method in [7] employs 158607 node points.

Next, we consider some PhC devices proposed in [7]. For the square lattice of dielectric rods used above, a microcavity can be formed by removing one single rod. In Fig. 6 (left panel), we show a microcavity coupled to two semi-infinite PhC waveguides. For a given incoming propagating mode in the left waveguide, we calculate the reflected and transmitted waves in the left and right waveguides, respectively. For this problem, we use a computation domain of 5×11 unit cells. In the vertical direction, we truncate the bulk PhC by retaining five rows of cylinders in each side of the defects, and use simple zero boundary conditions at the top and the bottom of the truncated domain. In the horizontal direction, the computation domain is bounded by the two dashed lines shown in Fig. 6. The boundary conditions on these two dashed lines are similar to (16) and (15), where *y* is now replaced by *x*. The computation domain involves 106 edges on which the field is to be determined. Therefore, the total number of unknowns is 106*N*, where *N* is the number of points on each edge of the unit cells. For *N*=5, we obtain the transmission spectrum of this structure as shown in Fig. 6 (right panel). We also consider a double microcavity coupled with waveguides, as shown in Fig. 7 (left panel). For this structure, our computation domain covers 8×11 unit cells and involves 179 edges. In Fig. 7 (right panel), we show the transmission spectrum obtained with *N*=5. For both single and double microcavities, our results agree fairly well with those reported in [7], except that the transmission spectra in [7] exhibit some oscillations. Resonant transmissions are observed at ω*a*/(2*πc*)=0.38672 for the single microcavity, and atωa/(2π*c*)=0.38415 and 0.38945 for the double microcavity.

Koshiba *et al* proposed some simple PhC waveguide branches in [7]. For the *Y*-branch shown in Fig. 8 (left panel), our computation domain contains 3×17 unit cells and 116 edges. Similar to the cases of microcavities coupled with waveguides, we use simple zero conditions at the top and bottom boundaries and rigorous boundary conditions, similar to (16) and (15), on the two vertical dashed lines. Using *N*=5 points on each edge of the unit cells, we obtain the transmission and reflection spectra shown in Fig. 8 (right panel), where the incoming wave is a propagating mode in the left waveguide (port 1). The *T*-branch shown in Fig. 9 is especially easy to analyze. Our computation domain involves only 17 unit cells and 50 edges. For *N*=5, we only have to solve a linear system for 250 unknowns and the results are shown in Fig. 9 (right panel). For both *Y*- and *T*-branches, our results and those reported in [7] are indistinguishable. Notice that the *T*-branch has a high transmission at frequencies satisfying 0.386≤ω*a*/(2*πc*)≤0.403.

In a recent paper [8], Ogusu and Takayama proposed PhC waveguides with stubs for possible applications as optical filters. In Fig. 10, we show the PhC waveguides with one or two stubs. The bulk PhC in Fig. 10 is a square lattice of dielectric rods in the air, where the refractive index and the radius of the rods are $\sqrt{11.9}$ and *r*=0.2*a*, respectively. The main waveguide is created by removing a row of rods. Each stub consists of two dielectric rods with a radius *r*
* _{s}*. For the waveguide with two stubs, the three rods between the stubs and closest to the waveguide core are also modified to have a radius

*r*

*. These structures were analyzed by the FDTD method in [8] using a computation domain of 11×17=187 unit cells surrounded by perfectly matched layers. In our calculations, the computation domains cover 11 and 55 unit cells for the cases with one and two stubs, respectively. The corresponding number of edges is 32 and 116. In Fig. 11, we show the transmission spectra of the waveguide with one stub for a few different values of*

_{e}*r*

*. For the waveguide with two stubs, we choose*

_{s}*r*

*=0.45*

_{s}*r*and consider a few different values of

*r*

*. The results are shown in Fig. 12. For both cases, our results and those reported in [8] agree in the main features, but differ in many details. It is likely that the computation domain used in their FDTD calculations is much too small. This is reflected in their transmission curve for single stub with*

_{e}*r*

*=*

_{s}*r*, as it deviates significantly from the constant 1. Our computation domains are even smaller, but the semi-infinite waveguides are rigorously simulated. The results in Fig. 11 and Fig. 12 are obtained with

*N*=7, therefore we only have to solve linear systems with 224 and 812 unknowns for waveguides with single and double stubs, respectively. To gain more confidence in our solutions, we have further increased the computation domain and the integer

*N*, but the results remain the same.

Finally, we consider an example previously analyzed by Ikuno and Naka [6]. As shown in Fig. 13 (left panel), the structure is a 60° waveguide bend with a microcavity, where the bulk PhC is a triangular lattice of dielectric rods in air. The waveguide is formed by removing one row of rods. The refractive index and the radius of the rods are 3.4 and 0.175*a*, respectively, where *a* is the lattice constant. The dielectric rod at the center of the microcavity is allowed to have a different refractive index. We analyze the transmission and reflection characteristics of the bend, assuming that the background PhC and the waveguides all extend to infinity and an incoming propagating mode is given in the horizontal waveguide. For this structure, we use hexagon unit cells and the computation domain shown in Fig. 13. The computation domain is obtained by truncating the bulk PhC in directions transverse to the waveguide axes with simple zero Dirichlet boundary conditions and terminating the semi-infinite waveguides with rigorous boundary conditions. Retaining five rows in each side of the waveguides and keeping some distance between the two boundaries for terminating the waveguides, we obtain the computation domain with 127 unit cells. For this structure, there are three distinct hexagon unit cells: the regular unit cell with a rod, the empty defect cell and the special cell at the center of the microcavity. For each unit cell, we calculate a (6*N*)×(6*N*) matrix approximation of the DtN map using *N* points on each edge. As described in section 4, the DtN maps of the unit cells are used to build the rigorous boundary conditions for terminating the waveguides. Finally, we set up a linear system of equations for the wave field on the cell edges. The total number of edges included in the linear system is 378 (about three edges per unit cell). In Fig. 13, we show the transmission spectra for a few different values of *ε*
* _{s}* (the dielectric constant of the special rod at the center) and for the

*E*polarization. Our results are obtained with

*N*=4 and they agree well with the FDTD results reported in [6]. For

*ε*

*=1, a resonant transmission is observed at ω*

_{s}*a*/(2

*πc*)=0.4084. In Fig. 14, we show the magnitude of the electric field at the resonant frequency. The field is particularly strong around the microcavity and it has nearly identical patterns in the two waveguides. Our method is efficient, since the linear system involves only 378×

*N*=1512 unknowns.

## 6. Conclusions

In this paper, we presented a Dirichlet-to-Neumann (DtN) map method for analyzing general two-dimensional photonic crystal (PhC) devices. The bulk PhC is a square or triangular lattice of circular cylinders in a homogeneous background medium. The device is connected by a few PhC waveguides (as ports) that extend to infinity. The method relies on the DtN maps of both regular and defect unit cells. A computation domain is obtained by terminating the PhC waveguides using rigorous boundary conditions. The problem is reduced to a relatively small linear system of equations on the edges of the unit cells in the computation domain. Accurate solutions are obtained using only 5 or 4 points on each edge for square or hexagon unit cells, respectively. This implies that we use only 10 or 12 unknowns for each unit cell in the final linear system of equations. In fact, these values correspond to representing the wave field in each unit cell by 20 or 24 cylindrical waves. For unit cells that are smaller than the free space wavelength, this is usually quite sufficient.

Most authors use the FDTD method to analyze PhC devices, but the accuracy may be limited if infinite PhC waveguides are truncated by perfectly matched layers. In earlier works, the DtN-map methodwas developed for eigenvalue problems associated with band structures [21,22,29], PhC waveguides [23] and microcavities [24], and to boundary values problems for finite PhCs between two homogeneous media [25–28]. The DtN-map method presented in this paper has a much wider scope as demonstrated in many numerical examples. Since the equation for an interior edge is only related to the edges of two neighboring unit cells (Section 2), the DtN-map method gives rise to sparse linear systems. In contrast, the multipole method [30–32] produces dense linear systems for coefficients in cylindrical wave expansions around the cylinders. However, the DtN-map method is currently limited to ideal two-dimensional structures that are invariant in the third dimension. Further work is needed to extend this method to devices in PhC slabs.

## Acknowledgments

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

## References and links

**1. **E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett. **58**, 2059–2062 (1987). [CrossRef] [PubMed]

**2. **S. John, “Strong localization of photons in certain disordered dielectric superlattices,” Phys. Rev. Lett. **58**, 2486–2489 (1987). [CrossRef] [PubMed]

**3. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn, *Photonic Crystals: Molding the Flow of Light*, (Princeton University Press, Princeton, NJ.1995).

**4. **A. Mekis, J. C. Chen, I. Kurland, S. H. Fan, P. R. Villeneuve, and J. D. Joannopoulos, “High transmission through sharp bends in photonic crystal waveguides,” Phys. Rev. Lett. **77**, 3787–3790 (1996). [CrossRef] [PubMed]

**5. **J. Smajic, C. Hafner, and D. Erni, “Design and optimization of an achromatic photonic crystal bend,” Opt. Express **11**, 1378–1384 (2003). [CrossRef] [PubMed]

**6. **H. Ikuno and Y. Naka, “Finite-difference time-domain method applied to photonic crystals,” in *Electromagnetic Theory and Appications for Photonic Crystals* , ed., K. Yasumoto, (CRC Press, Taylor & Francis Group, 2006).

**7. **M. Koshiba, Y. Tsuji, and M. Hikari, “Time-domain beam propagation method and its application to photonic crystal circuits,” J. Lightw. Technol. **18**, 102–110 (2000). [CrossRef]

**8. **K. Ogusu and K. Takayama, “Transmission characteristics of photonic crystal waveguides with stubs and their application to optical filters,” Opt. Lett. **32**, 2185–2187 (2007). [CrossRef] [PubMed]

**9. **T. Fujisawa and M. Koshiba, “Finite-element modeling of nonlinear interferometers based on photonic-crystal waveguides for all-optical signal processing,” J. Lightw. Technol. **24**, 617–623 (1006). [CrossRef]

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

**11. **L. C. Botten, N. A. Nicorovici, R. C. McPhedran, C. M. de Sterke, and A. A. Asatryan, “Photonic band structure calculations using scattering matrices,” Phys. Rev. E **64**, 046603 (2001). [CrossRef]

**12. **M. Marrone, V. F. Rodriguez-Esquerre, and H. E. Hernández-Figueroa, “Novel numerical method for the analysis of 2D photonic crystals: the cell method,” Opt. Express **10**, 1299–1304 (2002). [PubMed]

**13. **E. Moreno, D. Erni, and C. Hafner, “Band structure computations of metallic photonic crystals with the multiple multipole method,” Phys. Rev. B **65**, 155120 (2002). [CrossRef]

**14. **S. Jun, Y. S. Cho, and S. Im, “Moving least-square method for the band-structure calculation of 2D photonic crystals,” Opt. Express **11**, 541–551 (2003). [CrossRef] [PubMed]

**15. **C. P. Yu and H. C. Chang, “Compact finite-difference frequency-domain method for the analysis of two-dimensional photonic crystals,” Opt. Express **12**, 1397–1408 (2004). [CrossRef] [PubMed]

**16. **S. Guo, F. Wu, S. Albin, and R. S. Rogowski, “Photonic band gap analysis using finite-difference frequency-domain method,” Opt. Express **12**, 1741–1746 (2004). [CrossRef] [PubMed]

**17. **S. Y. Shi, C. H. Chen, and D. W. Prather, “Revised plane wave method for dispersive material and its application to band structure calculations of photonic crystal slabs,” Appl. Phys. Lett. **86**, 043104 (2005). [CrossRef]

**18. **S. Wilcox, L. C. Botten, R. C. McPhedran, C. G. Poulton, and C. M. de Sterke, “Modeling of defect modes in photonic crystals using the fictitious source superposition method,” Phys. Rev. E **71**, 056606 (2005). [CrossRef]

**19. **M. C. Lin and R. F. Jao, “Finite element analysis of photon density of states for two-dimensional photonic crystals with in-plane light propagation,” Opt. Express **15**, 207–218 (2007). [CrossRef] [PubMed]

**20. **P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E **75**, 026703 (2007). [CrossRef]

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

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

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

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

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

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

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

**30. **D. Felbacq, G. Tayeb, and D. Maystre, “Scattering by a random set of parallel cylinders,” J. Opt. Soc. Am. A **11**, 2526–2538 (1994). [CrossRef]

**31. **J. Yonekura, M. Ikeda, and T. Baba, “Analysis of finite 2-D photonic crystals of columns and lightwave devices using the scattering matrix method,” J. Lightw. Technol. **17**, 1500–1508 (1999). [CrossRef]

**32. **P. A. Martin, *Multiple Scattering*, (Cambridge University Press, Cambridge, UK, 2006). [CrossRef]