## Abstract

The Green’s functions are physical responses due to a single point source in a periodic lattice. The single point source can also correspond to an impurity or a defect. In this paper, the Green’s functions, including the scatterers, for periodic structures such as in photonic crystals and metamaterials are calculated. The Green’s functions are in terms of the multiband solutions of the periodic structures. The Green’s functions are broadband solutions so that the frequency or wavelength dependences of the physical responses can be calculated readily. They are obtained by integrating the periodic Green’s function including the scatterers in the Brillouin zone. Low wavenumber extraction methods are used to accelerate the convergence rate of the multiband expansions. The low wavenumber component represents the reactive near field. The band solutions of the periodic structure are obtained from a surface integral equation solution, which is effectively converted to a linear eigenvalue problem, giving multiple band solutions simultaneously. Numerical results are illustrated for the band modal fields, the periodic Green’s functions, and the single point source Green’s functions for two-dimensional (2D) perfect-electric-conductor (PEC) scatterers in a 2D lattice.

© 2017 Optical Society of America

## 1. INTRODUCTION

Waves in periodic structures are important in physics and engineering and in the design of photonic, electronic, acoustic, microwave, and millimeter devices such as those in photonic crystals, phononic crystals, and metamaterials [1–5]. In such problems of wave functional materials, there exists a lattice with scatterers embedded periodically in the lattice. The design of the lattice periodicity and the inclusions in a unit cell creates unique band structures with new wave phenomena such as edge states and topological insulators [6–8]. Besides the band structure, the band field solutions and the Green’s functions are also important physical quantities. The Green’s functions are physical responses of a point source. However, in periodic structures, only the periodic Green’s functions of an empty lattice have been studied extensively. In this paper, we calculate the Green’s functions of periodic structures including the scatterers, which can be of arbitrary shape [9–14]. The Green’s function offers physical insights into the wave behavior in the passband and stop band of photonic crystals and metamaterials. The Green’s function represents the response of sources and also impurities and defects. It can be used to formulate surface integral equations to deal with finite sizes, defects, imperfections, and impurities in the periodic lattice. The periodic Green’s function in an integral equation formulation has been used to derive the effective material parameters of the periodic structures [15–17]. The Green’s function including periodic scatterers has been constructed to study the dipole near field and radiation field inside a bandgap material [9–12]. The field excited by a dipole near a periodic structure has also been examined [13,14] using the periodic Green’s function.

The methods of band calculations include the plane-wave method [18–26] and the Korringa Kohn Rostoker (KKR) [27,28] or multiple-scattering (MS) method [29,30]. The plane-wave method casts the problem into a linear eigenvalue problem and provides multiple band solutions simultaneously for a point in the Brillouin zone. However, in the plane-wave method, the permittivity (or potential) profile is expanded in a Fourier series, making it difficult to use when there are abrupt changes in permittivity such as with scatterers, which are perfect-electric-conductors (PECs) [21]. For dielectric scatters, smoothed and effective permittivities are introduced near the dielectric interfaces to improve the convergence [22,23,26]. The KKR method formulates surface integral equations with the periodic Green’s function of an empty lattice. For spherical objects, the scattering T-matrix can be calculated analytically at a single wavenumber, and the scattering can be included in the boundary condition [29,30]. The method is only convenient for scatterers with shapes of separable geometries such as circular cylinders or spheres. In the KKR/MS method, the eigenvalue problem is nonlinear, requiring the root-seeking procedure of calculating the modal wavenumber, one root at a time. For scatterers of arbitrary shape, numerical approaches such as the finite-difference time-domain method [31,32] have been used.

Using the concept of modal expansion of the periodic Green’s function, we have recently developed the method of broadband Green’s function with low wavenumber extraction (BBGFL) [33–39] that gives an accelerated convergence of the multiple band expansions. Using BBGFL, surface integral equations are formulated and solved by the method of moment (MoM) so that the method is applicable to scatterers of arbitrary shape. The determination of modal band solutions in this method is a linear eigenvalue problem, so the multiband solutions are computed for a Bloch wavenumber simultaneously. This is in contrast to using the usual free-space Green’s function or the KKR/MS method in which the eigenvalue problem is nonlinear. The modal field solutions are wavenumber independent. We have applied the BBGFL to calculate band diagrams of periodic structures. The BBGFL method is applicable to both PEC [33] and dielectric periodic scatterers [34]. The method is broadband so that the frequency or wavelength dependences can be calculated readily. Our method has some similarities to the hybrid-plane-wave- and integral-equation-based methods [40–42], where an integral–differential eigensystem is derived for an auxiliary extended problem that has smooth eigenfunctions.

The goal of this paper is to calculate the Green’s function for periodic structures that includes infinite periodic scatterers and to illustrate physically the responses due to point sources in the periodic structures. Such Green’s functions also correspond to responses due to an impurity or defect in the periodic structures. The mathematical steps are (1) to solve for the band modal fields and normalize the band modal fields [35,37], and (2) to calculate the periodic Green’s function at a single low wavenumber ${k}_{L}$ from surface integral equation and (3) the periodic Green’s function at any wavenumber $k$ using the accelerated modal representation for each Bloch wave-vector in the first Brillouin zone, and (4) to calculate the Green’s function due to a single point source by integrating the periodic Green’s function over the first Brillouin zone with respect to the Bloch wavenumber [9–14]. Our approach is related to [9–12] by representing the periodic Green’s function in terms of multiband solutions and in applying the phased-array method to obtain the point source response. However, we apply the surface integral equation with the MoM to solve for the multiple band solutions instead of using the plane-wave expansion, making the approach applicable to high permittivity contrast, arbitrary shape scatterers, and nonpenetrable scatterers. We use the low wavenumber extraction technique to accelerate the convergence of the band modal representation, making a broadband response ready to obtain. The low wavenumber component represents the reactive near field, which is the origin of slow convergence.

The outline of the paper is as follows. In Section 2, we formulate the problem and state several Green’s functions. In Section 3, we derive the broadband Green’s functions for a single point source. In Section 4, we illustrate numerical results of the fields on the surfaces of the scatterers, fields of band solutions, and the Green’s functions including scatterers. The band solutions and the calculation of the periodic Green’s function at a single low wavenumber are described in Appendices A and B, respectively.

## 2. BROADBAND GREEN’S FUNCTION WITH LOW WAVENUMBER EXTRACTION BASED ON MODAL REPRESENTATION

Consider a 2D periodic lattice, as illustrated in Fig. 1, in the $x\u2013y$ plane with $z$ perpendicular to the plane. We use 2D vector $\overline{\rho}=x\widehat{x}+y\widehat{y}$. The lattice vector is ${\overline{R}}_{pq}=p{\overline{a}}_{1}+q{\overline{a}}_{2}$, where $-\infty <p,q<\infty $. The lattice has periodic scatterers as in photonic crystals and metamaterials. We use superscript “0” to denote the Green’s function for an empty lattice and superscript “$S$” to denote the Green’s function and the band modal field with the scatterers. We use ${\overline{k}}_{i}$ to denote a wave vector in the first Brillouin zone. Let ${\mathrm{\Omega}}_{0}=|{\overline{a}}_{1}\times {\overline{a}}_{2}|$ be the size of the unit cell. The reciprocal lattice vectors are ${\overline{K}}_{mn}=m{\overline{b}}_{1}+n{\overline{b}}_{2}$, where ${\overline{b}}_{1}=2\pi \frac{{\overline{a}}_{2}\times \widehat{z}}{{\mathrm{\Omega}}_{0}}$, and ${\overline{b}}_{2}=2\pi \frac{\widehat{z}\times {\overline{a}}_{1}}{{\mathrm{\Omega}}_{0}}$. The point ${\overline{k}}_{i}$ in the first Brillouin zone is ${\overline{k}}_{i}({\beta}_{1},{\beta}_{2})={\beta}_{1}{\overline{b}}_{1}+{\beta}_{2}{\overline{b}}_{2}$, $0\le {\beta}_{1},{\beta}_{2}\le 1$.

The periodic Green’s function is the response to an infinite array of periodic point sources with progressive phase shift, denoted by ${\delta}^{\infty}({\overline{k}}_{i};\overline{\rho}-{\overline{\rho}}^{\prime})$, as given by

The broadband Green’s function is based on band modal representation. It is broadband because the modal field distributions are independent of wavenumber. The broadband Green’s function has a simple dependence on wavenumber $k$, so the Green’s functions can be calculated at any $k$ readily. The low wavenumber extraction at a single wavenumber ${k}_{L}$ accelerates the convergence of the band modal summation. With the extraction, the band modal summation converges when the field point coincides with the source point. The low wavenumber extraction requires a single MoM solution at that ${k}_{L}$. There are three BBGFLs in this paper.

- (a) The usual empty lattice periodic Green’s function ${g}_{P}^{0}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ at ${\overline{k}}_{i}$ of the first Brillouin zone but with low wavenumber extraction. The empty lattice ${g}_{P}^{0}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ is due to periodic sources of Eq. (1).
- (b) The periodic Green’s function including periodic scatters ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ at ${\overline{k}}_{i}$ of the first Brillouin zone. The ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ is also due to periodic sources of Eq. (1).
- (c) The Green’s function including the periodic scatterers ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime})$. ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime})$ is due to a single source at ${\overline{\rho}}^{\prime}$, that is $\delta (\overline{\rho}-{\overline{\rho}}^{\prime})$. Thus, ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime})$ is nonperiodic.

In this section, we outline the steps of how each of the three BBGFLs is calculated. More details are given in subsequent sections, the appendices, and [37].

- (a) The empty lattice Green’s function ${g}_{P}^{0}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$.

In the spectral domain, it is

This can be viewed as modal expansion, where $|{\overline{k}}_{i\alpha}|$ is the modal frequency and the Bloch wave ${\tilde{\psi}}_{\alpha}^{0}({\overline{k}}_{i},\overline{\rho})=\frac{1}{\sqrt{{\mathrm{\Omega}}_{0}}}\text{\hspace{0.17em}}\mathrm{exp}({\overline{k}}_{i\alpha}\xb7\overline{\rho})$ is the normalized modal field. The modal field satisfies the orthonormal condition.

The above expansion is of poor convergence and does not converge at $\overline{\rho}={\overline{\rho}}^{\prime}$. In previous papers [33,34], we use a low wavenumber extraction by first calculating ${g}_{p}^{0}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ as

- (b) The periodic Green’s function including the periodic scatterers ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$.

Using BBGFL ${g}_{P}^{0}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$, we formulated a surface integral equation and solved it with the MoM to find the band modal solutions. Using BBGFL in the MoM, we obtain a linear eigenvalue problem to calculate simultaneously the band solutions with all the band eigenvalues ${k}_{\beta}^{S}({\overline{k}}_{i})$ dependent on ${\overline{k}}_{i}$ and not on $k$. In Appendix A, we describe efficient methods of calculating band modal solution ${\psi}_{\beta}^{S}({\overline{k}}_{i},\overline{\rho})$ and perform the normalization of modal fields, which is

We next use the normalized modal solutions ${\psi}_{\beta}^{S}({\overline{k}}_{i},\overline{\rho})$ to calculate the Green’s function ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$, which includes the presence of the scatterers. Note that ${k}_{\beta}^{S}({\overline{k}}_{i})$ is real. We have

Equation (5) has poor convergence. Thus, we also use a single low wavenumber extraction ${k}_{L}$ to obtain

Note that ${g}_{P,B}^{S}(k,{k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ converges with respect to $1/{[{k}_{\beta}^{S}({\overline{k}}_{i})]}^{4}$ in contrast to $1/{[{k}_{\beta}^{S}({\overline{k}}_{i})]}^{2}$ as in Eq. (5). Only relatively few modes are needed to construct the broadband Green’s function. In truncating the series of Eq. (6), we require $[{k}_{\beta}^{S}({\overline{k}}_{i})]\gg k$, i.e., the maximum $|{k}_{\beta}^{S}({\overline{k}}_{i})|$ included should be several times larger than the largest $k$ of interest to ensure convergence. In Appendix B, we describe the method of calculating ${g}_{P}^{S}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ at a single ${k}_{L}$.

- (c) The Green’s function including the periodic scatterers ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime})$ due to a single source at ${\overline{\rho}}^{\prime}$.

It will be shown in Section 3 that

It is noteworthy that enhanced convergence techniques [38,39] are under active development, ready to improve the convergence rate of the broadband Green’s function of Eqs. (3) and (6) from fourth power to sixth power, further reducing the number of modes required to accurately construct the Green’s function. This also enables the broadband Green’s function for applications in 3D vector wave equations [39].

## 3. GREEN’S FUNCTION ${\mathit{g}}^{\mathit{S}}(\mathit{k};\overline{\mathit{\rho}},{\overline{\mathit{\rho}}}^{\prime})$ OF A SINGLE POINT SOURCE: INTEGRATION OVER THE BRILLOUIN ZONE

We find the Green’s function ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime})$ due to a single point source $\delta (\overline{\rho}-{\overline{\rho}}^{\prime})$ in the lattice including periodic scatterers. This can be obtained by integration of the periodic Green’s function ${g}_{P}^{S}(k,{\overline{k}}_{i}({\beta}_{1},{\beta}_{2});\overline{\rho},{\overline{\rho}}^{\prime})$ of a periodic point source array of ${\delta}^{\infty}(\overline{\rho}-{\overline{\rho}}^{\prime};{\overline{k}}_{i})$, as given in Eq. (1). Note that ${\overline{k}}_{i}$ is a function of ${\beta}_{1}$ and ${\beta}_{2}$, as defined in Section 2. In this section, we derive the relations between ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime})$ and ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$.

It can be shown that $\delta (\overline{\rho}-{\overline{\rho}}^{\prime})$ can be represented by integrating ${\delta}^{\infty}({\overline{k}}_{i};\overline{\rho}-{\overline{\rho}}^{\prime})$ over the first Brillouin zone [9–14],

In showing the above, we make use of the property ${\overline{a}}_{i}\xb7{\overline{b}}_{j}=2\pi {\delta}_{ij}$. It then follows that

The integrand is a periodic function with respect to ${\overline{k}}_{i}$. We apply the following midpoint rectangular quadrature rule for its numerical evaluation [14]:

## 4. RESULTS AND DISCUSSION

#### A. Green’s Function ${\mathit{g}}_{\mathit{P}}^{\mathit{S}}(\mathit{k},{\overline{\mathit{k}}}_{\mathit{i}};\overline{\mathit{\rho}},{\overline{\mathit{\rho}}}^{\prime \prime})$ at a Single ${\mathit{k}}_{\mathit{L}}$

We illustrate results considering a periodic array of circular PEC cylinders. The primary lattice vectors are defined by ${\overline{a}}_{1}=\frac{a}{2}(\sqrt{3}\widehat{x}+\widehat{y})$, and ${\overline{a}}_{2}=\frac{a}{2}(-\sqrt{3}\widehat{x}+\widehat{y})$, where $a=1$ is the normalized lattice constant. The cylinders have radii of $b=0.2a$ centering at ${\overline{R}}_{pq}$, as defined in Section 2. The background region out of the cylinders has permittivity of ${\u03f5}_{b}=8.9{\u03f5}_{0}$. The source point is at ${\overline{\rho}}^{\prime \prime}=\frac{1}{3}({\overline{a}}_{1}+{\overline{a}}_{2})$. We choose ${k}_{L}=\frac{2\pi}{a}{f}_{nL}\sqrt{\frac{{\u03f5}_{b}}{{\u03f5}_{0}}}$, where ${f}_{nL}=0.001$. The normalized frequency ${f}_{n}$ is defined as ${f}_{n}=\frac{ka}{2\pi}\sqrt{\frac{{\u03f5}_{0}}{{\u03f5}_{b}}}$. The Brillouin zone vector ${\overline{k}}_{i}({\beta}_{1},{\beta}_{2})$ is as defined in Section 2. We choose ${\beta}_{1}=0.1$ and ${\beta}_{2}=0.05$. The position of the cylinder and the source point in the unit lattice are illustrated in Fig. 2(a). We also choose a special field point at $\overline{\rho}=\frac{7}{16}{\overline{a}}_{1}$, where we will examine the ${\overline{k}}_{i}$ dependence. The surface currents on the cylinder at ${f}_{nL}$ are shown in Fig. 2(b), which demonstrate a peak at $\varphi =90\xb0$, closest to the source point.

The field distribution of ${g}_{P}^{S}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ over the lattice is shown in Fig. 3. The field repeats itself under Bloch wave conditions. The repeating of the point sources is clear from the magnitude of the field distribution. The phase progression according to ${\overline{k}}_{i}$ is manifested in the real and imaginary parts of the field distributions.

#### B. Results of Band Modal Fields ${\mathit{\psi}}_{\mathit{\beta}}^{\mathit{S}}({\overline{\mathit{k}}}_{i};\overline{\mathit{\rho}})$

We use the same periodic array of PEC cylinders to illustrate the spatial variation of the modal fields. We also choose the same ${\overline{k}}_{i}$ and ${k}_{L}$. In Fig. 4, we show the modal fields for the lowest three modes. The field solutions are extinguished inside the PEC cylinder. They exhibit more complicated patterns over the lattice as the normalized modal frequency increases. The field patterns are orthogonal to each other.

#### C. Illustration of Results of ${\mathit{g}}_{\mathit{P}}^{\mathit{S}}(\mathit{k},{\overline{\mathit{k}}}_{\mathit{i}};\overline{\mathit{\rho}},{\overline{\mathit{\rho}}}^{\prime \prime})$ at a General Wavenumber $k$

In this subsection, we show the results of ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ at a general wavenumber $k$. We compute according to the broadband expression of Eq. (7) using the modal fields ${\psi}_{\beta}^{S}({\overline{k}}_{i};\overline{\rho})$ that are independent of $k$. The expression is applicable to all $k$. To test the correctness of the expression, we can also compare the results by computing directly using the approach of Appendix B for a general $k$.

We first illustrate the results of ${g}_{P,B}^{S}(k,{k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ following Eq. (6) at ${f}_{n}=0.2$ using all the modes with ${k}_{\beta}\le {k}_{\mathrm{max}}$. Equation (6) suggests a convergence rate of up to $1/{k}_{\beta}^{4}$ subject to ${k}_{\beta}\gg k$. The convergence is checked by taking ${k}_{\mathrm{max}}=8k$ as the benchmark and calculating the relative error in evaluating ${g}_{P,B}^{S}$ as the number of modes included increases. The test confirms the rapid convergence of ${g}_{P,B}^{S}$. In this case, truncating the series with ${k}_{\beta}/k\le 3$ yields a relative error of less than 1% ([37], Section 7.3.3). It is noted that the convergence rate in general slows down as $k$ increases, and more modes are required in the series to ensure accuracy. The density of modes also increases as the frequency increases. In the following computations, we set ${k}_{\mathrm{max}}=8k$, which includes 49 modes for ${f}_{n}=0.2$. The results of ${g}_{P,B}^{S}$ are given in Fig. 5. Comparing to the results of ${g}_{P}^{S}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$, as shown in Fig. 3, the behavior of ${g}_{P,B}^{S}(k,{k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ is smooth and without singularity, unlike ${g}_{P}^{S}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$. Figure 5 also shows that the spatial variation of ${g}_{P,B}^{S}$ is close to the modal field distribution of ${\tilde{\psi}}_{\beta}^{S}({\overline{k}}_{i};\overline{\rho})$ at ${f}_{n\beta}=0.216$, as shown in Fig. 4. This is because ${f}_{n}$ is close to ${f}_{n\beta}$ so that the contribution from the corresponding mode becomes dominant in shaping the broadband Green’s function ${g}_{P,B}^{S}(k,{k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ at this $k$.

We then use Eq. (7) to evaluate ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ at three different values of $k$ corresponding to ${f}_{n}=0.1$, 0.2, and 0.4, respectively. The spatial variations of ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ are plotted in Fig. 6. The modes included are with ${k}_{\beta}\le {k}_{\mathrm{max}}=8k$. The results show that the Green’s function at a specific ${\overline{k}}_{i}$ varies significantly with $k$.

We report the CPU times required to evaluate ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ following Eq. (7) in Table 1, where the absolute time unit is ignored. The CPU times are decomposed into several steps. The modal analysis takes up 10.5 time units to solve Eq. (A9) and obtain all the modal wavenumber ${k}_{\beta}^{S}$ and the eigenvectors ${\overline{b}}_{\beta}$. The eigenvector ${\overline{b}}_{\beta}$ is related to the modal field ${\psi}_{\beta}^{S}({\overline{k}}_{i};\overline{\rho})$ through Eq. (A20). The calculation of ${\psi}_{\beta}^{S}({\overline{k}}_{i};\overline{\rho})$ using Eq. (A20) on 100 field points $\overline{\rho}$ takes up 0.2425 time units. The calculation of ${g}_{P,B}^{S}(k,{k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ through Eq. (6) on 100 field points takes up 0.0905 time units. The evaluation of ${g}_{P}^{S}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ at the single low wavenumber uses up 6.748 time units to assemble the impedance matrix and to solve for the surface currents. This part of CPU can be saved almost completely if the same ${k}_{L}$ is chosen as in the first step of modal analysis. Another 0.5414 time units are used to calculate ${g}_{P}^{S}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ at 100 field points. Note that only the time fraction spent on calculating ${g}_{P,B}^{S}(k,{k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ duplicates when using Eq. (7) to evaluate ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ on multiple values of $k$. This does not introduce a considerable time increase as compared to using Eqs. (B1) and (B2) to evaluate ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ at each $k$. Using Eqs. (B1) and (B2) leads to linearly scaling CPU times with respect to the number of frequencies considered. The benefits of using Eq. (7) over Eqs. (B1) and (B2) become more obvious as the number of frequencies to be examined increases. The gain in efficiency approaches $\sim 80.5$ for a large enough number of frequencies.

In Fig. 7, we plot ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ as a function of the normalized frequency ${f}_{n}$ at field point $\overline{\rho}=\frac{7}{16}{\overline{a}}_{1}$ and source point ${\overline{\rho}}^{\prime \prime}=\frac{1}{3}({\overline{a}}_{1}+{\overline{a}}_{2})$. We show the results of two cases: (a) lossless background with ${\u03f5}_{b}=8.9{\u03f5}_{0}$ and (b) lossy background with ${\u03f5}_{b}=8.9(1+0.11i){\u03f5}_{0}$. In both cases, we have used a real low wavenumber ${k}_{L}$ corresponding to ${f}_{nL}=0.001$. The results of Eq. (7) with BBGFL are compared to the results of Appendix B by solving the surface integral equation (SIE) directly. For the lossless case, the results agree well, except close to ${f}_{n}=0.22$, which is close to the modal frequency of ${f}_{n\beta}=0.216$. The poles in the modal expansion of the Green’s function as exhibited in the factor of $1/({k}_{\beta}^{2}-{k}^{2})$ in Eqs. (5) and (6) causes the decrease in accuracy. The errors can be due to numerical errors in evaluating the modal wavenumber ${k}_{\beta}$, which exaggerates the disagreement as $k\to {k}_{\beta}$. The errors can be also due to errors in evaluating the modal fields ${\psi}_{\beta}^{S}$ and due to insufficient modes included in the modal expansion of Eq. (6). Improved convergence techniques [38,39] are being developed to further reduce the number of modes required in the expansion and improve the accuracy in evaluating ${k}_{\beta}$ and ${\psi}_{\beta}^{S}$. The SIE results also become less accurate near resonance. The agreement is improved in the lossy case since by using a complex $k$, the resonance issue is avoided. We also examine the differences from the SIE results. The differences are generally within 5%, except close to modal frequencies. The results of the lossy case are also more accurate.

#### D. Illustration of Results of the Green’s Function ${\mathit{g}}^{\mathit{S}}(\mathit{k};\overline{\mathit{\rho}},{\overline{\mathit{\rho}}}^{\prime \prime})$ of a Single Point Source

We now examine the performance of BBGFL when used to calculate ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime \prime})$ over a wide spectrum after integrating ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ over the Brillouin zone.

We first illustrate the integrand of Eq. (10) ${g}_{P}^{S}(k,{\overline{k}}_{i}({\beta}_{1},{\beta}_{2});\overline{\rho},{\overline{\rho}}^{\prime \prime})$ as a function of ${\overline{k}}_{i}$ for the stop band and the passband in Figs. 8 and 9, respectively. We also decompose the integrand of ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ into the primary contribution ${g}_{P}^{0}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ and the response contribution ${g}_{P}^{R}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$. The results are evaluated at ${f}_{n}=0.2$ in stop band, with source point ${\overline{\rho}}^{\prime \prime}=\frac{1}{3}({\overline{a}}_{1}+{\overline{a}}_{2})$ and field point $\overline{\rho}=\frac{7}{16}{\overline{a}}_{1}$. A lossless background with ${\u03f5}_{b}=8.9{\u03f5}_{0}$ is assumed. It is noted that in Figs. 8(b) and 8(c), both the primary and response components change rapidly as a function of ${\overline{k}}_{i}$. The results are that the singular parts cancel each other, such that ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ is a smooth function of ${\overline{k}}_{i}$, as shown in Fig. 8(a).

We note that ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ is smooth in Fig. 8(a) because the frequency is in the stop band so that there are no poles/modes encountered over the entire Brillouin zone. In Fig. 9, we show the passband results at ${f}_{n}=0.26$. In Fig. 9(a), we plot the integrand of ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ with a lossless background ${\u03f5}_{b}=8.9{\u03f5}_{0}$. The integrand is singular when the modes ${k}_{\beta}$ coincide with $k$. In Fig. 9(b), ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ is shown with a lossy background ${\u03f5}_{b}=8.9(1+0.11i){\u03f5}_{0}$. The complex $k$ avoids the real pole of ${k}_{\beta}$ giving smooth results.

In Fig. 10, we show the results of integration over the Brillouin zone and plot the spatial variations of ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime \prime})$ in a lossless background at ${f}_{n}=0.1$, ${f}_{n}=0.2$, and ${f}_{n}=0.4$, respectively. The results are calculated following Eq. (10) invoking the low wavenumber extraction technique. To simplify, we have chosen a constant ${f}_{nL}=0.001$ over the entire Brillouin zone. Comparing to the band diagram of the periodic structure given in Ref. [33], we see the spreading of the fields at ${f}_{n}=0.4$ as it is in the passband.

In Fig. 11, we plot ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime \prime})$ as a function of the normalized frequency. The results are evaluated at source point ${\overline{\rho}}^{\prime \prime}=\frac{1}{3}({\overline{a}}_{1}+{\overline{a}}_{2})$ and field point $\overline{\rho}=\frac{7}{16}{\overline{a}}_{1}$. The values of the Green’s function obtained with BBGFL are compared to the results by substituting in Eq. (10) the ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ obtained by solving the integral equations of Eqs. (B2) and (B1) for each $k$ and ${\overline{k}}_{i}$. In Fig. 11(a), we assume a lossless background of ${\u03f5}_{b}=8.9{\u03f5}_{0}$. Note that other than the peak value at ${f}_{n}=0.26$, the agreement is good. The oscillation of the Green’s function is closely related to the band diagrams of the periodic structure. It is suppressed in the stop band below ${f}_{n}=0.2$ and becomes complex above that frequency. One should keep in mind the results based on the direct SIE are also subject to errors when entering the passband due to the poles in ${g}_{P}^{S}$. In Fig. 11(b), we assume a lossy background permittivity of ${\u03f5}_{b}=8.9(1+0.11i){\u03f5}_{0}$. This yields a complex $k$, and avoids the poles in the modal expansion. The agreement for the lossy case is better.

The spatial variations of ${g}^{S}(k;\overline{\rho},{\overline{\rho}}^{\prime \prime})$ with the lossy background are plotted in Fig. 12 at the three normalized frequencies of ${f}_{n}=0.1$, ${f}_{n}=0.2$, and ${f}_{n}=0.4$, respectively. Comparing with Fig. 10, the spread-out of field in the passband of ${f}_{n}=0.4$ is reduced due to absorption loss.

## 5. CONCLUSIONS

In this paper, we calculate and illustrate the Green’s functions due to point sources inside a periodic array of scatterers. The forms of the Green’s function are in terms of band solutions and are broadband so that they can be evaluated over wide spectrum of frequencies or wavelengths. The Green’s functions provide a physical understanding of the propagation and scattering in periodic structures. We have illustrated Green’s functions in the bandgap and in the passband. We are presently using this Green’s function to formulate integral equations that can be used to model excitations, impurities, displacement of scatterers, disorder, defects, and finite-sized periodic structures [39]. Extensions to the 3D case are also presently studied [39].

## APPENDIX A: BAND ANALYSIS USING BBGFL—${\mathit{k}}_{\mathit{\beta}}^{\mathit{S}}({\overline{\mathit{k}}}_{\mathit{i}})$ AND ${\mathit{\psi}}_{\mathit{\beta}}^{\mathit{S}}({\overline{\mathit{k}}}_{i};\overline{\mathit{\rho}})$

## 1. Calculation of ${\mathit{k}}_{\mathit{\beta}}^{\mathit{S}}({\overline{\mathit{k}}}_{\mathit{i}})$ and ${\mathit{\psi}}_{\mathit{\beta}}^{\mathit{S}}({\overline{\mathit{k}}}_{\mathit{i}};\overline{\mathit{\rho}})$ Using BBGFL

We summarize the procedures that we developed in Refs. [33,34] to characterize the modal wavenumber ${k}_{\beta}^{S}({\overline{k}}_{i})$ and modal field ${\psi}_{\beta}^{S}({\overline{k}}_{i};\overline{\rho})$ using BBGFL. We follow a convention of using subscript $\beta $ to denote the modes including the periodic scatterers, while using $\alpha $ to denote the Bloch modes of the empty lattice. We limit ourselves to the case of PEC scatterers with TMz polarization.

The surface integral equation to solve is

Applying the MoM with pulse basis and point matching, Eq. (A2) is converted into

where ${q}_{n}=\mathrm{\Delta}{t}_{n}J({\overline{\rho}}_{n})$, $n=\mathrm{1,2},\dots ,N$, $\mathrm{\Delta}{t}_{n}$ and ${\overline{\rho}}_{n}$ are the arc length and the center of each small patch, respectively, and $N$ is the number of discretization in representing the boundary. $\overline{\overline{R}}$ is a matrix of $N\times M$, where $M$ is the number of Bloch waves included in the broadband Green’s function, which also determines the dimension of the eigenvalue problem to be solved. ${R}_{m\alpha}={R}_{\alpha}({\overline{\rho}}_{m})$. $\overline{\overline{L}}$ is a matrix of $N\times N$, where ${L}_{mn}$ is the impedance matrix elements evaluated at ${k}_{L}$. We haveAfter discretization, Eq. (A3) can be also converted into a matrix form as

We can represent $\overline{q}$ in terms of $\overline{b}$ from Eq. (A5) as

and then use Eq. (A8) to eliminate $\overline{q}$ in Eq. (A7) and convert it into the following linear eigenvalue problem: where $\overline{b}$ appear as the eigenvectors, and $\lambda =1/({k}^{2}-{k}_{L}^{2})$ is the eigenvalue corresponding to modal wavenumber ${k}_{\beta}^{S}({\overline{k}}_{i})$. The matrix $\overline{\overline{A}}=\overline{\overline{D}}-{\overline{\overline{R}}}^{\u2020}{\overline{\overline{L}}}^{-1}\overline{\overline{R}}$ is of dimenstion $M\times M$.The modal field ${\psi}_{\beta}^{S}({\overline{k}}_{i};\overline{\rho})$ is then calculated through the following extinction theorem:

By using BBGFL, the solution of the linear eigenvalue problem Eq. (A9) yields all the eigenmodes and modal fields simultaneously.

## 2. Normalization of Band Modal Solutions

For the modal expansion of Eq. (5) to hold, the modal field ${\psi}_{\beta}^{S}({\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$ must satisfy the orthonormal condition of Eq. (4). In this section, we describe the efficient procedure of making the modal solutions orthonormal [35].

Substituting Eq. (3) into Eq. (A10),

Using the wave equations that govern ${\psi}_{\beta}^{S}(\overline{\rho};{\overline{k}}_{i})$, ${\tilde{\psi}}_{\alpha}^{0}({\overline{k}}_{i};\overline{\rho})$, and ${g}_{P}^{0}({k}_{L},{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime})$,

Substitute Eq. (A11) into the right-hand side of Eq. (A12) and making use of Eqs. (A13) and (A14),

Using the condition of ${k}_{L}\to 0$, i.e., ${k}_{L}\ll {k}_{\beta}^{S}$ and ${k}_{L}\ll |{\overline{k}}_{i\alpha}|$, where both ${k}_{\beta}^{S}$ and ${\overline{k}}_{i\alpha}$ depend on ${\overline{k}}_{i}$, Eq. (A15) simplifies to

The normalized modal field is then

After Eq. (A16) and invoking the orthonormal condition of ${\tilde{\psi}}_{\alpha}^{0}(\overline{\rho};{\overline{k}}_{i})$, it follows that

We assume that $\sum _{\alpha}{b}_{\alpha \beta}^{*}{b}_{\alpha {\beta}^{\prime}}={\delta}_{\beta {\beta}^{\prime}}$ following the Sturm–Liouville theory. Also, they can be orthonormalized through a Gram–Schmidt process. Thus,

Thus, the normalization is

We will then use ${\psi}_{\beta}^{S}$ to denote ${\tilde{\psi}}_{\beta}^{S}$ and assume that it is normalized.

The normalized modal field can be calculated from either Eqs. (A11) and (A17) or Eq. (A20). Note that Eq. (A20) is only valid for ${k}_{L}\to 0$.

## APPENDIX B: GREEN’S FUNCTION ${\mathit{g}}_{\mathit{P}}^{\mathit{S}}(\mathit{k},{\overline{\mathit{k}}}_{\mathit{i}};\overline{\mathit{\rho}},{\overline{\mathit{\rho}}}^{\prime \prime})$ AT A SINGLE LOW WAVENUMBER ${\mathit{k}}_{\mathit{L}}$

We formulate the surface integral equations to calculate ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$. The surface integral equation is formulated for a general $k$. However, in implementation, we only need to solve for a single low wavenumber ${k}_{L}$.

Let ${\overline{\rho}}^{\prime \prime}$ be the source point outside the scatterer and $\overline{\rho}$ be the field point outside the scatterer, which can also approach the surface from outside the scatterer; see Fig. 1. The extinction theorem governing ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ is then given by

We use the Dirichlet boundary condition ${g}_{P}^{S}(k,{\overline{k}}_{i};{\overline{\rho}}^{\prime},{\overline{\rho}}^{\prime \prime})=0$. This is the case when we examine the electric field response due to a $\widehat{z}$-polarized line source outside of a periodic array of PEC scatterers. The fields are TM polarized to $z$.

By letting $\overline{\rho}$ approach the boundary ${S}_{00}$, we obtain the following SIE:

The SIE in Eq. (B2) is to solve the surface current $J({\overline{k}}_{i};{\overline{\rho}}^{\prime})={\widehat{n}}^{\prime}\xb7{\nabla}^{\prime}{g}_{P}^{S}(k,{\overline{k}}_{i};{\overline{\rho}}^{\prime},{\overline{\rho}}^{\prime \prime})$ at a single low wavenumber $k={k}_{L}$.

Notice that this is the same equation that governs the modal analysis problem as we developed in Refs. [33,34]. The only difference is replacing the right-hand side (excitation) with the direct incidence field from the periodic point source array ${g}_{P}^{0}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$. We apply the pulse basis and point matching in the MoM. The evaluation of the matrix elements and the right-hand side also follow directly the scheme that we developed in Refs. [33,34]. Note that solving the surface integral equation can have multiple right-hand sides by varying ${\overline{\rho}}^{\prime \prime}$ as the MoM impedance matrix remains the same.

After the surface currents are obtained, we use Eq. (B1) to obtain the low wavenumber Green’s function ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ with $k={k}_{L}$.

Note that both ${g}_{P}^{0}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ and ${g}_{P}^{S}(k,{\overline{k}}_{i};\overline{\rho},{\overline{\rho}}^{\prime \prime})$ satisfy the Bloch wave condition

## REFERENCES

**1. **J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, *Photonic Crystals: Molding the Flow of Light* (Princeton University, 2011).

**2. **S. Tretyakov, *Analytical Modeling in Applied Electromagnetics* (Artech House, 2003).

**3. **J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, “Magnetism from conductors and enhanced nonlinear phenomena,” IEEE Trans. Microwave Theory Tech. **47**, 2075–2084 (1999). [CrossRef]

**4. **D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, “Composite medium with simultaneously negative permeability and permittivity,” Phys. Rev. Lett. **84**, 4184–4187 (2000). [CrossRef]

**5. **S. Yang, J. H. Page, Z. Liu, M. L. Cowan, C. T. Chan, and P. Sheng, “Focusing of sound in a 3D phononic crystal,” Phys. Rev. Lett. **93**, 024301 (2004). [CrossRef]

**6. **A. B. Khanikaev, S. H. Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald, and G. Shvets, “Photonic topological insulators,” Nat. Mater. **12**, 233–239 (2013). [CrossRef]

**7. **Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong, and B. Zhang, “Topological acoustics,” Phys. Rev. Lett. **114**, 114301 (2015). [CrossRef]

**8. **Z. Wang, Y. D. Chong, J. D. Joannopoulos, and M. Soljacic, “Reflection-free one-way edge modes in a gyromagnetic photonic crystal,” Phys. Rev. Lett. **100**, 013905 (2008). [CrossRef]

**9. **T. Suzuki and P. K. I. Yu, “Emission power of an electric dipole in the photonic band structure of the FCC lattice,” J. Opt. Soc. Am. B **12**, 570–582 (1995). [CrossRef]

**10. **C. Caloz, D. Curcio, A. Alvarez-Melcon, A. K. Skrivervik, and F. E. Gardiol, “Slot antenna on a photonic crystal substrate Green’s function study,” Proc. SPIE **3795**, 176–187 (1999). [CrossRef]

**11. **C. Caloz, A. K. Skrivervik, and F. E. Gardiol, “An efficient method to determine Green’s functions of a two-dimensional photonic crystal excited by a line source—the phased-array method,” IEEE Trans. Microwave Theory Tech. **50**, 1380–1391 (2002). [CrossRef]

**12. **M. G. M. V. Silveirinha and C. A. Fernandes, “Radiation from a short vertical dipole in a disk-type PBG material,” in *IEEE Antennas and Propagation Society International Symposium* (IEEE, 2003), Vol. 3, pp. 990–993.

**13. **F. Capolino, D. R. Jackson, and D. R. Wilton, “Fundamental properties of the field at the interface between air and a periodic artificial material excited by a line source,” IEEE Trans. Antennas Propag. **53**, 91–99 (2005). [CrossRef]

**14. **F. Capolino, D. R. Jackson, D. R. Wilton, and L. B. Felsen, “Comparison of methods for calculating the field excited by a dipole near a 2-D periodic material,” IEEE Trans. Antennas Propag. **55**, 1644–1655 (2007). [CrossRef]

**15. **M. Silveirinha and C. Fernandes, “Effective permittivity of metallic crystals: a periodic Green’s function formulation,” Electromagnetics **23**, 647–663 (2003). [CrossRef]

**16. **M. G. Silveirinha, “Metamaterial homogenization approach with application to the characterization of microstructured composites with negative parameters,” Phys. Rev. B **75**, 115104 (2007). [CrossRef]

**17. **M. G. Silveirinha, “Generalized Lorentz-Lorenz formulas for microstructured materials,” Phys. Rev. B **76**, 245117 (2007). [CrossRef]

**18. **K. M. Ho, C. T. Chan, and C. M. Soukoulis, “Existence of a photonic gap in periodic dielectric structures,” Phys. Rev. Lett. **65**, 3152–3155 (1990). [CrossRef]

**19. **K. M. Leung and Y. F. Liu, “Full vector wave calculation of photonic band structures in face-centered-cubic dielectric media,” Phys. Rev. Lett. **65**, 2646–2649 (1990). [CrossRef]

**20. **R. D. Mead, K. D. Brommer, A. M. Rappe, and J. D. Joannopoulos, “Existence of a photonic bandgap in two dimensions,” Appl. Phys. Lett. **61**, 495–497 (1992). [CrossRef]

**21. **H. S. Sozuer, J. W. Haus, and R. Inguva, “Photonic bands: convergence problems with the plane-wave method,” Phys. Rev. B **45**, 13962–13972 (1992). [CrossRef]

**22. **R. D. Meade, A. M. Rappe, K. D. Brommer, J. D. Joannopoulos, and O. L. Alerhand, “Accurate theoretical analysis of photonic band-gap materials,” Phys. Rev. B **48**, 8434–8437 (1993). [CrossRef]

**23. **S. G. Johnson, “Erratum: accurate theoretical analysis of photonic band-gap materials,” Phys. Rev. B **55**, 15942 (1997).

**24. **M. Plihal and A. A. Maradudin, “Photonic band structure of two-dimensional systems: the triangular lattice,” Phys. Rev. B **44**, 8565–8571 (1991). [CrossRef]

**25. **V. Kuzmiak, A. A. Maradudin, and F. Pincemin, “Photonic band structures of two-dimensional systems containing metallic components,” Phys. Rev. B **50**, 16835–16844 (1994). [CrossRef]

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

**27. **J. Korringa, “On the calculation of the energy of a Bloch wave in a metal,” Physica **13**, 392–400 (1947). [CrossRef]

**28. **W. Kohn and N. Rostoker, “Solution of the Schrödinger equation in periodic lattices with an application to metallic lithium,” Phys. Rev. **94**, 1111–1120 (1954). [CrossRef]

**29. **K. M. Leung and Y. Qiu, “Multiple-scattering calculation of the two-dimensional photonic band structure,” Phys. Rev. B **48**, 7767–7771 (1993). [CrossRef]

**30. **W. Zhang, C. T. Chan, and P. Sheng, “Multiple scattering theory and its application to photonic band gap systems consisting of coated spheres,” Opt. Express **8**, 203–208 (2001). [CrossRef]

**31. **S. Fan, P. R. Villeneuve, and J. D. Joannopoulos, “Large omnidirectional band gaps in metallodielectric photonic crystals,” Phys. Rev. B **54**, 11245–11251 (1996). [CrossRef]

**32. **A. J. Ward and J. B. Pendry, “Calculating photonic Green’s function using a nonorthogonal finite-difference time-domain method,” Phys. Rev. B **58**, 7252–7259 (1998). [CrossRef]

**33. **L. Tsang, “Broadband calculations of band diagrams in periodic structures using the broadband Green’s function with low wavenumber extraction (BBGFL),” Prog. Electromagn. Res. **153**, 57–68 (2015). [CrossRef]

**34. **L. Tsang and S. Tan, “Calculations of band diagrams and low frequency dispersion relations of 2D periodic dielectric scatterers using broadband Green’s function with low wavenumber extraction (BBGFL),” Opt. Express **24**, 945–965 (2016). [CrossRef]

**35. **L. Tsang and S. Huang, “Broadband Green’s function with low wave number extraction for arbitrary shaped waveguide with applications to modeling of vias in finite power/ground plane,” Prog. Electromagn. Res. **152**, 105–125 (2015). [CrossRef]

**36. **S. Huang and L. Tsang, “Fast electromagnetic analysis of emissions from printed circuit board using broadband Green’s function method,” IEEE Trans. Electromagn. Compat. **58**, 1642–1652 (2016). [CrossRef]

**37. **S. Tan, “Multiple volume scattering in random media and periodic structures with applications in microwave remote sensing and wave functional materials,” Ph.D. dissertation (University of Michigan, 2016).

**38. **T.-H. Liao, K.-H. Ding, and L. Tsang, “High order extractions of broadband Green’s function with low wavenumber extractions for arbitrary shaped waveguide,” Prog. Electromagn. Res. **158**, 7–20 (2017). [CrossRef]

**39. **S. Tan and L. Tsang, “Scattering from finite periodic arrays using broadband Green’s function of periodic scatterers with low wavenumber extraction (BBGFL),” in IEEE Progress in Electromagnetics Research Symposium (PIERS), St. Petersburg, Russia, May 22–25, 2017.

**40. **M. G. Silveirinha and C. A. Fernandes, “Efficient calculation of the band structure of artificial materials with cylindrical metallic inclusions,” IEEE Trans. Microwave Theory Tech. **51**, 1460–1466 (2003). [CrossRef]

**41. **M. G. Silveirinha and C. A. Fernandes, “A hybrid method for the efficient calculation of the band structure of 3-D metallic crystals,” IEEE Trans. Microwave Theory Tech. **52**, 889–902 (2004). [CrossRef]

**42. **M. G. Silverinha and C. A. Fernandes, “Computation of the electromagnetic modes in two-dimensional photonic crystals: a technique to improve the convergence rate of the plane wave method,” Electromagnetics **26**, 175–187 (2006). [CrossRef]