Accurate analysis of electromagnetic scattering from periodic circular cylinder array with defects

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 transitionmatrix 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 Optical Society of America OCIS codes: (050.0050) Diffraction and gratings; (050.1755) Computational electromagnetic methods; (050.1950) Diffraction gratings. 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). 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). 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). 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). 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). 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). 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). 11. K. Watanabe, J. Pištora, and Y. Nakatake, “Coordinate transformation formulation of electromagnetic scattering from imperfectly periodic surfaces,” Opt. Express (to be published). #162163 $15.00 USD Received 26 Jan 2012; revised 13 Mar 2012; accepted 23 Mar 2012; published 24 Apr 2012 (C) 2012 OSA 7 May 2012 / Vol. 20, No. 10 / OPTICS EXPRESS 10646 12. N. A. Nicorovici and R. C. McPhedran, “Lattice sums for off-axis electromagnetic scattering by gratings,” Phys. Rev. E 50, 3143–3160 (1994). 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). 14. C. A. Balanis, Advanced Engineering Electromagnetics (Wiley, New York, 1989).


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

Settings of the problem
The present paper considers the scattering problem of electromagnetic fields with a timedependence 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 μ c , and their radius is denoted by 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 ε s and the permeability μ s .
The wavenumber and the characteristic impedance in each medium are respectively denoted by k j = ω (ε j μ j ) 1/2 and ζ j = (μ j /ε j ) 1/2 for j = c, s.The axis of the lth-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 D, which is a finite subset of the integer set Z. If an integer l is an element of D, the cylinder whose center is at (x, y) = (l d, 0) is removed.Also, the complement of D in Z is denoted by D c .The electromagnetic fields are supposed to be uniform in the 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.

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 g g (Z) (x, y), in which the nth-component is given as with where 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 q ) and (x r , y r ) be the reference points of the bases.Then, when (x, y) is inside a circle with center (x r , y r ) and radius ρ(x q − x r , y q − y r ), the translation relation for the cylindrical-wave bases is given by where G G 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 where the superscript t denotes the transpose and a a a q 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 where a a a (s) q denotes the column matrix generated by the expansion coefficients of the scattered waves from the qth-cylinder.Using the translation relation Eq. ( 4), the total field near but outside the qth-cylinder (q ∈ D c ) can be expressed as where D c \{q} denotes the set in which an element q is removed from D c .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 qth-cylinder.Since the relation between the coefficient matrices of the incident and the scattered fields are given by the transitionmatrix (T-matrix), we have a a a (s) The (n, m)-entries of the T-matirx T T T are given by for TM-polarization for TE-polarization (10) where the symbol δ n,m stands for Kronecker's delta.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 {a a a (s) q } 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 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 with The column matrices a a a (i) (ξ ) and a a a (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 d .The original coefficient matrices are inversely obtained by integrating on the transform parameter ξ as a a a for f = i, s.Substituting Eq. ( 9) into Eq.( 16) and using Eqs.( 15) and ( 17), we obtain the following relation: with where the superscript −1 stands for the matrix inverse.The entries of matrix L L 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.
Since the transformed field ψ (s,p) (x; ξ , y) can be expressed in the same form with Eq. ( 14), the coefficient matrix a a a (s) (ξ ) is also decomposed as a a a (s) (ξ ) = a a a (s,p) (ξ ) + a a a (s,d) where a a a (s,p) (ξ ) and a a a (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 (D = / 0), the matrices M M 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 Substituting Eqs. ( 23) and (24) into Eq.( 18), we may obtain the following equation for a a a (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 {ξ l } L l=1 in the Brillouin zone −k d /2 < ξ ≤ k d , 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 l=1 .Then, we may derive the following relation: with From Eq. ( 26), the coefficient matrices of the residual field at the sample point is obtained by

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 Using the Fourier integral representation for the Hankel function of the first kind of order zero [5], the nth-components of the column matrices a a a (i) (ξ ) are given by a a a (i) (ξ ) where the infinite summation in Eq. ( 31) is also truncated in practical computation.We use 2 K + 1 terms (from −Kthto Kth-order terms) to obtain the following results.Also, we use the same sample points {ξ l } L l=1 for a a a (i) (ξ ), and the column matrix b b b l given in Eq. ( 28) are then approximated as 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 (a) (C) 2012 OSA 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 q ) = ψ(x p , y p ; x q , y q ) − ψ(x q , y q ; x p , u p ) ψ(x p , y p ; x q , y q ) (35) where ψ(x p , y p ; x q , y q ) denotes the field observed at (x p , y p ) for a line-source located at (x q , y q ).The reciprocity theorem requires that this function is zero when both (x p , y p ) and (x q , y q ) are 0 Line-Source located in the surrounding medium.We fix one point at (x p , y p ) = (0, 2 d) and the other point (x q , y q ) is moved on the lines 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.

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 with k x = −k s cos(θ (i) ) and k y = −k s sin(θ (i) ).Then, the column matrices a a a (i) (ξ ), a a a (s,p) (ξ ), and b b b l are given by a a a (i) where I I I denotes the identity matrix and the nth-components of the column matrix a a a (i) 0 are given by a a a 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.

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.(a) 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 (C) 2012 OSA 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.

FieldFig. 2 .
Fig. 2. Convergence of the field intensities at (x, y) = (0, ±d) for a line-source excitation.The dotted curves are obtained by the trapezoidal scheme, and the solid curves are obtained by the Gauss-Legendre scheme applying to split subintervals.(a) Convergence on the number of sample points L (b) Convergence on the truncation order K.

FieldFig. 6 .Fig. 7 .Fig. 8 .
Fig. 6.Convergence of the field intensities at (x, y) = (0, ±d) for the plane-wave incidence.The dotted curves are obtained by the trapezoidal scheme, and the solid curves are obtained by the Gauss-Legendre scheme applying to split subintervals.(a) Convergence on the number of sample points L (b) Convergence on the truncation order K.