## Abstract

This paper considers the two-dimensional electromagnetic scattering from periodic array of circular cylinders in which some cylinders are removed, and presents a formulation based on the recursive transition-matrix algorithm (RTMA). The RTMA was originally developed as an accurate approach to the scattering problem of a finite number of cylinders, and an approach to the problem of periodic cylinder array was then developed with the help of the lattice sums technique. This paper introduces the concept of the pseudo-periodic Fourier transform to the RTMA with the lattice sums technique, and proposes a spectral-domain approach to the problem of periodic cylinder array with defects.

© 2012 OSA

## 1. Introduction

The photonic crystals have been proposed as media in which photons can be manipulated. These structures can exhibit a complete photonic bandgap. A lot of complete-photonic-bandgap based applications, functionality is provided through the precise, controlled incorporation of three-dimensionally engineered defects [1, 2]. The photonic crystals with multiple-hole defects have been compared with a single hole defect photonic crystal structures to demonstrate the advantage of using the multiple-hole defects for sensing applications [3]. In laser technology, the decimated photonic crystal lattices play important role [4]. These lattices have reduced number of holes, which enhance the current and heat conductive paths from the optical cavity. Many of these photonic crystals consist of arrays of dielectric circular cylinders, and more effective technique is often required to analyze the electromagnetic fields in the cylinder array.

The recursive transition-matrix algorithm (RTMA) [5, 6] is one of the most commonly used approaches to the electromagnetic scattering from a finite number of parallel cylinders. It uses the cylindrical-wave expansions to express the field components and the boundary conditions at the cylinder surfaces are derived using Graf’s addition theorem. Roussel et al. [7] extended the RTMA to the plane-wave scattering from periodic cylinder array, which consists of infinite number of cylinders. When the plane-wave illuminates a periodic structure, the Floquet theorem asserts that the scattered fields are pseudo-periodic (namely, each field component is given by a product of a periodic function and an exponential phase factor, and has an equidistant discrete spectrum in the wavenumber space). Many of the numerical approaches to periodic structures use this property, and the formulation proposed by Roussel et al. is also based on the Floquet theorem. If the incident field is not the plane-wave or the structure is not perfectly periodic, the fields generally have continuous spectra in the wavenumber space, and an artificial discretization in the wavenumber space is necessary for the spectral-domain approach. Watanabe and Yasumoto introduced the pseudo-periodic Fourier transform (PPFT) [8] to consider the discretization scheme for periodic structures for non-plane incident waves, and this approach has applied to the scattering problem of the circular cylinders located near the periodic cylinder array [9]. This paper considers the electromagnetic scattering problem of the periodic circular cylinder array in which some cylinders are removed, and presents a formulation based on the RTMA with the help of PPFT concept. Since the PPFT makes the fields pseudo-periodic, the conventional grating theories with the use of the generalized Fourier series expansions also become possible to be applied to the scattering problems of imperfectly periodic structures. The combinations with the rigorous coupled-wave analysis and with the coordinate transformation method were recently proposed in separated papers [10, 11].

## 2. Settings of the problem

The present paper considers the scattering problem of electromagnetic fields with a time-dependence exp(−*iωt*) from a periodic circular cylinder array with defects. Figure 1 shows an example of the structures under consideration. The structure consists of identical circular cylinders, which are infinitely long in the *z*-direction and situated parallel to each other. The cylinders are made with a homogeneous and isotropic medium described by the permittivity *ɛ _{c}* and the permeability

*μ*, and their radius is denoted by

_{c}*a*. For a lossy medium, we use complex values for the permittivity and/or the permeability, and the terms concerning to the current densities are eliminated in the following formulation. The surrounding region is filled by a lossless, homogeneous, and isotropic medium with the permittivity

*ɛ*and the permeability

_{s}*μ*. The wavenumber and the characteristic impedance in each medium are respectively denoted by

_{s}*k*=

_{j}*ω*(

*ɛ*)

_{j}μ_{j}^{1/2}and

*ζ*= (

_{j}*μ*/

_{j}*ɛ*)

_{j}^{1/2}for

*j*=

*c,s*. The axis of the

*l*th-cylinder is located at (

*x,y*) = (

*l d,*0) for integer

*l*, though some cylinders are removed from the periodic array. To indicate the removed cylinders, we introduce a notation

*𝒟*, which is a finite subset of the integer set . If an integer

*l*′ is an element of

*𝒟*, the cylinder whose center is at (

*x,y*) = (

*l*′

*d,*0) is removed. Also, the complement of

*𝒟*in is denoted by

*𝒟*. The electromagnetic fields are supposed to be uniform in the

^{c}*z*-direction, and two-dimensional scattering problem is here considered. Two fundamental polarizations are expressed by TM and TE, in which the magnetic and the electric fields are respectively perpendicular to the

*z*-axis. We denote the

*z*-component of the electric field for the TM-polarization and the

*z*-component of the magnetic field for the TE-polarization by

*ψ*(

*x,y*), and show the formulation for both polarizations simultaneously. The incident field

*ψ*

^{(i)}(

*x,y*) is supposed to illuminate the cylinders from outside and there exists no source inside the cylinders.

## 3. Formulation

The present formulation uses mainly the cylindrical-wave expansions to express the fields in the surrounding medium. The basis functions for the cylindrical-wave expansion are here given by column matrices *g*^{(Z)}(*x,y*), in which the *n*th-component is given as

*Z*specifies the cylindrical functions associating with the cylindrical-wave bases in such a way that

*Z*=

*J*denotes the Bessel function and

*Z*=

*H*

^{(1)}denotes the Hankel function of the first kind. Let (

*x*,

_{q}*y*) and (

_{q}*x*,

_{r}*y*) be the reference points of the bases. Then, when (

_{r}*x, y*) is inside a circle with center (

*x*,

_{r}*y*) and radius

_{r}*ρ*(

*x*), the translation relation for the cylindrical-wave bases is given by

_{q}− x_{r}, y_{q}− y_{r}

*G*^{(Z)}(

*x,y*) denotes the Toeplitz matrix whose (

*n,m*)-entries are given by

We have supposed that there exists no sources inside the inhomogeneous region −*a* ≤ *y* ≤ *a* Therefore, if we choose the reference point at (*x,y*) = (*qd*,0) for an integer *q*, the incident field in the surrounding medium can be expanded in a series of the cylindrical-waves associated with the Bessel functions as

*t*denotes the transpose and ${\mathit{a}}_{q}^{(i)}$ is a column matrix generated by the expansion coefficients. On the other hand, the scattered field

*ψ*

^{(s)}(

*x,y*) outside the cylinders consists of the outward propagating waves from the cylinders, and it can be expressed as

*q*th-cylinder. Using the translation relation Eq. (4), the total field near but outside the

*q*th-cylinder (

*q*∈

*𝒟*) can be expressed as

^{c}*𝒟*\{

^{c}*q*} denotes the set in which an element

*q*is removed from

*𝒟*. The first and the second terms on the right-hand side of Eq. (8) are given by superpositions of cylindrical-waves associated with the Bessel function and the Hankel function of the first kind, respectively, and they represent the incident and the scattered fields for the

^{c}*q*th-cylinder. Since the relation between the coefficient matrices of the incident and the scattered fields are given by the transition-matrix (T-matrix), we have

*n,m*)-entries of the T-matirx

**are given by**

*T**δ*stands for Kronecker’s delta.

_{n,m}The process up to here is same with that of the conventional RTMA [5, 6]. The conventional method considers the problem of a finite number of cylinders and the coefficient matrices {
${\mathit{a}}_{q}^{(s)}$} are numerically obtained from Eq. (9) by truncating the expansions. However, the structure under consideration includes an infinite number of cylinders, and Eq. (9) cannot be directly solved. To resolve this difficulty, we introduce the PPFT. Let *f*(*x*) be a function of *x*. Then the PPFT and its inverse transform are formally defined by

*ξ*is the transform parameter and

*k*= 2

_{d}*π/d*denotes the inverse lattice constant. The transformed functions have a pseudo-periodic property in terms of

*x: f̄*(

*x − qd*;

*ξ*) =

*f̄*(

*x*;

*ξ*) exp(−

*iqdξ*) for any integer

*q*, and also have a periodic property in terms of

*ξ*:

*f̄*(

*x*;

*ξ*

*− qk*) =

_{d}*f̄*(

*x*;

*ξ*). Applying the PPFT to the incident and the scattered fields, the transformed fields are expressed as follows

*ā*^{(i)}(

*ξ*) and

*ā*^{(s)}(

*ξ*) give the coefficients of the cylindrical-wave expansions for the reference point at (

*x,y*) = (0,0), and they are periodic in terms of

*ξ*with the period

*k*. The original coefficient matrices are inversely obtained by integrating on the transform parameter

_{d}*ξ*as

*f*=

*i,s*.

Substituting Eq. (9) into Eq. (16) and using Eqs. (15) and (17), we obtain the following relation:

**(**

*L**ξ*) are called the lattice sums [12] and known to converge very slowly. Yasumoto and Yoshitomi [13] have proposed an integral representation of the lattice sums, which makes it possible to obtain the accurate values of lattice sums in practical computation. By contrast, the scalar function

*c*

^{(d)}(

*ξ, ξ*′) is given by a sum of finite terms, and it expresses the effects of defects in the periodic array.

The scattered field *ψ*^{(s)}(*x,y*) is here decomposed into the scattered field by the periodic cylinder array without defects *ψ*^{(s,p)}(*x,y*) and the residual field *ψ*^{(s,d)}(*x,y*) = *ψ*^{(s)}(*x,y*) − *ψ*^{(s,p)}(*x,y*). Since the transformed field *ψ̄*^{(s,p)}(*x*;*ξ*,*y*) can be expressed in the same form with Eq. (14), the coefficient matrix *ā*^{(s)}(*ξ*) is also decomposed as

*ā*^{(s,p)}(

*ξ*) and

*ā*^{(s,d)}(

*ξ*) denote the coefficient matrices of

*ψ̄*

^{(s,p)}(

*x*;

*ξ*,

*y*) and

*ψ̄*

^{(s,d)}(

*x*;

*ξ,y*), respectively. If there is no defect in the cylinder array ( $\mathcal{D}=\overline{)0}$), the matrices

*M*^{(d)}(

*ξ*,

*ξ*′) and

*c*

^{(d)}(

*ξ,ξ*′) vanish. Therefore, from Eq. (18), the coefficient matrix for the scattered field by the periodic cylinder array is given as

*ā*^{(s,d)}(

*ξ*):

To solve the integral Eq. (25), we introduce here a discretization in the transform parameter *ξ*. Considering the periodicity, we take *L* sample points
${\{{\xi}_{l}\}}_{l=1}^{L}$ in the Brillouin zone −*k _{d}*/2 <

*ξ*≤

*k*, and assume that Eq. (25) is satisfied at these points. Also, the integral on the left-hand side is approximated by an appropriate numerical integration scheme with the use of same sample points. We denote the respective weights of the sample points by ${\{{w}_{l}\}}_{l=1}^{L}$. Then, we may derive the following relation:

_{d}## 4. Numerical experiments

To validate the present formulation, we consider a specific example of periodic circular cylinder array with defects. The parameters are chosen as follows: *d* = 0.8 *λ*_{0}, *a* = 0.4 *d*, *ɛ _{c}* = 4

*ɛ*

_{0},

*ɛ*=

_{s}*ɛ*

_{0},

*μ*=

_{c}*μ*=

_{s}*μ*

_{0}, and three cylinders with the center at (

*x, y*) = (−2

*d*, 0), (0, 0), (2

*d,*0) are removed from the perfectly periodic array (

*𝒟*= {−2, 0, 2}). When implementing a practical computation, the cylindrical-wave expansions must be truncated. We denote the truncation order for the expansions by

*K*that truncates the expansions from −

*K*th- to

*K*th-order.

#### 4.1. Line-source excitation

First, we consider the fields excited by a line-source of unit amplitude, which is situated parallel to the *z*-axis at (*x,y*) = (*x*_{0}, *y*_{0}) for *y*_{0} > *a* or *y*_{0} < −*a*. The incident field is then expressed as

*n*th-components of the column matrices

*ā*^{(i)}(

*ξ*) are given by

*K*+ 1 terms (from −

*K*th- to

*K*th-order terms) to obtain the following results. Also, we use the same sample points ${\{{\xi}_{l}\}}_{l=1}^{L}$ for

*ā*^{(i)}(

*ξ*), and the column matrix

*b**given in Eq. (28) are then approximated as*

_{l}We set the location of the line-source at (*x*_{0}, *y*_{0}) = (*d*, 2*d*), and the obtained field intensities at (*x,y*) = (0, ±*d*) are shown in Fig. 2. Figure 2(a) shows the convergence characteristics in terms of the number of sample points *L*. The truncation order is set to *K* = 4 for calculating these results. The dotted curves are the results of the trapezoidal scheme, which uses equidistant sample points and equal weights. They oscillate and converge very slowly. The solid curves are obtained by the discretization scheme used in Refs. [8, 9]. In this scheme, the Brillouin zone is split at the Wood-Rayleigh anomalies and the Gauss-Legendre scheme is applied for each subinterval to decide the sample points and weights. Figure 2(a) shows that this scheme is effective to improve the convergence though the convergence for the TM-polarization is some slower than that for the TE-polarization. Figure 2(a) shows the convergence in terms of the truncation order for cylindrical-wave expansion *K*. The values are computed with *L* = 80, and the sample points and the weights are determined by applying the Gauss-Legendre scheme for the subintervals. The convergence with respect to the truncation order *K* is very fast. From the physical point of view, the results of the present formulation are not expected to noticeably different in many cases from those of the conventional RTMA [5] for the scattering from a finite number of cylinders, if the cylinder number is large enough and the observation points are located near the line-source. Here, we consider 98 cylinders located at (*x,y*) = (*qd*,0) for *q* = ±1, ±3, ±4, ±5,...,±50, and the field intensities at (*x,y*) = (0,*d*) are calculated by the conventional RTMA. The obtained values are 0.338 for the TM-polarization and 0.264 for the TE-polarization. These values are in good agreement with the results shown in Fig. 2. The total field intensities near the defects are computed with *L* = 80 and *K* = 4 by changing the observation point, and shown in Fig. 3. The positions of cylinder surfaces are indicated by the white dashed lines and the obtained results seem to be proper. Using the cylindrical coordinate (*ρ*^{(o)}, *ϕ*^{(o)}) and applying the saddle-point method for *ρ*^{(o)} → ∞ [14], an expression of the scattered field in the far-zone is obtained, and Fig. 4 shows the absolute values of the scattering pattern function.

We examine also the reciprocal property by using the reciprocity error defined by

*ψ*(

*x*,

_{p}*y*;

_{p}*x*,

_{q}*y*) denotes the field observed at (

_{q}*x*,

_{p}*y*) for a line-source located at (

_{p}*x*,

_{q}*y*). The reciprocity theorem requires that this function is zero when both (

_{q}*x*,

_{p}*y*) and (

_{p}*x*,

_{q}*y*) are located in the surrounding medium. We fix one point at (

_{q}*x*,

_{p}*y*) = (0, 2

_{p}*d*) and the other point (

*x*,

_{q}*y*) is moved on the lines

_{q}*y*= ±

*d*. The reciprocity errors at 101 equidistant points in −8

*d*≤

*x*≤ 8

*d*are calculated with

*L*= 80 and

*N*= 4 in the standard double-precision arithmetic, and the obtained results are shown in Fig. 5. The obtained values are smaller than 3 × 10

^{−14}and, considering the computation precision, it can be said that the reciprocity relation is perfectly satisfied.

#### 4.2. Plane-wave incidence

Next, we consider the plane-wave incidence problem. The angle of incidence *θ*^{(i)} is measured counter-clockwise from the *x*-axis, and the incident field of unit amplitude is expressed as

*k*= −

_{x}*k*cos(

_{s}*θ*

^{(i)}) and

*k*= −

_{y}*k*sin(

_{s}*θ*

^{(i)}). Then, the column matrices

*ā*^{(i)}(

*ξ*),

*ā*^{(s,p)}(

*ξ*), and

*b**are given by*

_{l}**denotes the identity matrix and the**

*I**n*th-components of the column matrix ${\mathit{a}}_{0}^{(i)}$ are given by

Figure 6 shows the results of the same convergence tests with Fig. 2 but for the angle of incidence *θ*^{(i)} = 70°. The convergence characteristics are similar to those for the line-source excitation, though the amplitudes of the oscillations appeared on the dotted curves are smaller. The results of the present formulation are compared with the those of the conventional RTMA for 98 cylinders located at (*x,y*) = (*qd,*0) for *q* = ±1, ±3, ±4, ±5,...,±50. The values obtained by the conventional RTMA are 1.377 for the TM-polarization and 0.417 for the TE-polarization. These values are in good agreement with the results shown in Fig. 6. The total field intensities near the defects and the scattering pattern of the residual field *ψ*^{(s,d)} are computed with *L* = 80 and *K* = 4, and shown in Figs. 7 and 8, respectively. In this computation, the Brillouin zone is split at the Wood-Rayleigh anomalies, and the sample points and the weights are determined by applying the Gauss-Legendre scheme for the subintervals.

## 5. Concluding remarks

The paper has presented an accurate method for analyzing the two-dimensional electromagnetic scattering from a periodic circular cylinder array in which some cylinders are removed. The present formulation was based on the RTMA with the lattice sums technique and the PPFT. We have shown numerical results for the line-source excitation and the plane-wave incidence. A simple discretization scheme presented in Ref. [8, 9] showed good convergences, and the numerical results are validated by comparing with that for an array of a large number of cylinders. The present formulation was also validated by the reciprocity test. The numerical results shown in this paper were for the array of lossless dielectric circular cylinders. However, the convergence characteristics for the array of lossy cylinders are similar in many cases. Also, the arrays of the perfectly conducting cylinders or the non-circular cylinders are possible to be treated by replacing the T-matrix, because only the T-matrix depends on the material and the shape of cylinders in the present formulation. When the cylinders are made of an anisotropic material or the fields are not uniform in the direction of cylinder axes, the TM- and the TE-polarizations cannot be decomposed, but the extension is straightforward.

## Acknowledgments

This work has been partially supported by the Grant Agency of the Czech Republic (#P205/11/2137), by Structural Funds of the European Union and state budget of the Czech Republic (NANOBASE project #CZ.1.07/2.3.00/20.0074), and by the IT4Innovations Centre of Excellence project, reg. no. CZ.1.05/1.1.00/02.0070.

## References and links

**1. **S. A. Rinne, F. García-Santamaría, and P. V. Braun, “Embedded cavities and waveguides in three-dimensional silicon photonic crystals,” Nat. Photonics **2**, 52–56 (2007). [CrossRef]

**2. **J. Ouellette, “Seeing the future in photonic crystals,” Ind. Phys . **7**, 14–17 (2001).

**3. **Ch. Kang and S. M. Weiss, “Photonic crystal with multi-hole defect for sensor applications,” Opt. Express **16**, 18188–18193 (2008). [CrossRef] [PubMed]

**4. **A. V. Giannopoulos, J. D. Sulkin, Ch. M. Long, J. J. Coleman, and K. D. Choquette, “Decimated photonic crystal defect cavity lasers,” IEEE J. Sel. Top. Quantum Electron . **17**, 1693–1694 (2011). [CrossRef]

**5. **W. C. Chew, *Waves and Fields in Inhomogeneous Media* (Van Nostrand Reinhold, New York, 1990).

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

**7. **H. Roussel, W. C. Chew, F. Jouvie, and W. Tabbara, “Electromagnetic scattering from dielectric and magnetic gratings of fibers — a T-matrix solution,” J. Electromagn. Waves Appl . **10**, 109–127 (1996). [CrossRef]

**8. **K. Watanabe and K. Yasumoto, “Two-dimensional electromagnetic scattering of non-plane incident waves by periodic structures,” Prog. Electromagn. Res . **PIER 74**, 241–271 (2007). [CrossRef]

**9. **K. Watanabe and Y. Nakatake, “Spectral-domain formulation of electromagnetic scattering from circular cylinders located near periodic cylinder array,” Prog. Electromagn. Res. B **31**, 219–237 (2011).

**10. **K. Watanabe, J. Pištora, and Y. Nakatake, “Rigorous coupled-wave analysis of electromagnetic scattering from lamellar grating with defects,” Opt. Express **19**, 25799–25811 (2011). [CrossRef]

**11. **K. Watanabe, J. Pištora, and Y. Nakatake, “Coordinate transformation formulation of electromagnetic scattering from imperfectly periodic surfaces,” Opt. Express(to be published). [PubMed]

**12. **N. A. Nicorovici and R. C. McPhedran, “Lattice sums for off-axis electromagnetic scattering by gratings,” Phys. Rev. E **50**, 3143–3160 (1994). [CrossRef]

**13. **K. Yasumoto and K. Yoshitomi, “Efficient calculation of lattice sums for free-space periodic Green’s function,” IEEE Trans. Antennas Propag . **47**, 1050–1055 (1999). [CrossRef]

**14. **C. A. Balanis, *Advanced Engineering Electromagnetics* (Wiley, New York, 1989).