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–y$ 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$.

 

Fig. 1. Illustration of periodic scatterers in 2D periodic lattice in the $x–y$ plane. $S_{pq}$ denotes the surface of the $pq$th scatterer. ${\overline{\rho }}^{\prime \prime }$ and $\overline{\rho }$ represent the locations of arbitrary source and field points, respectively.

Download Full Size | PDF

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.

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

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}}\exp ({\overline{k}}_{i\alpha }·\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

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

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·{\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 ${ϵ}_b=8.9{ϵ}_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{{ϵ}_b}{{ϵ}_0}}$, where $f_{nL}=0.001$. The normalized frequency $f_n$ is defined as $f_n=\frac{ka}{2\pi }\sqrt{\frac{{ϵ}_0}{{ϵ}_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°$, closest to the source point.

 

Fig. 2. (a) Geometry of the cylinder (red circle) and the source point (black cross) inside the unit cell. Blue circles denote $16\times 16$ uniformly distributed grid points where we probe the fields. The black plus sign denotes a special field point to be examined more closely. (b) Magnitude of surface currents on the cylinder.

Download Full Size | PDF

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.

 

Fig. 3. Field distribution of $g_P^S(k_L,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ over the lattice: (a) magnitude, (b) real part, and (c) imaginary part.

Download Full Size | PDF

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.

 

Fig. 4. Modal field distribution (magnitude) for the lowest three modes: (a) $f_n=0.216$, (b) $f_n=0.368$, and (c) $f_n=0.413$.

Download Full Size | PDF

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_{\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_{\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_{\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$.

 

Fig. 5. Spatial variation of $g_{P,B}^S(k,k_L,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ at $f_n=0.2$ with $f_{nL}=0.001$ and ${\overline{k}}_i=0.1{\overline{b}}_1+0.05{\overline{b}}_2$: (a) magnitude, (b) real part, and (c) imaginary part.

Download Full Size | PDF

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_{\max }=8k$. The results show that the Green’s function at a specific ${\overline{k}}_i$ varies significantly with $k$.

 

Fig. 6. Magnitude of $g_P^S(k,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ at (a) $f_n=0.1$, (b) $f_n=$ 0.2, and (c) $f_n=$ 0.4. The numbers of modes included in $g_B^S$ are 12, 49, and 116, respectively.

Download Full Size | PDF

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.

Table 1. CPU Time Decomposition

View Table

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 ${ϵ}_b=8.9{ϵ}_0$ and (b) lossy background with ${ϵ}_b=8.9(1+0.11i){ϵ}_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.

 

Fig. 7. $g_P^S(k,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ as a function of the normalized frequency $f_n$: (a) ${ϵ}_b=8.9{ϵ}_0$ and (b) ${ϵ}_b=8.9(1+0.11i){ϵ}_0$.

Download Full Size | PDF

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 ${ϵ}_b=8.9{ϵ}_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).

 

Fig. 8. Magnitude of the integrand as a function of ${\overline{k}}_i$: (a) $g_P^S(k,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$, (b) $g_P^0(k,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$, and (c) $g_P^R(k,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$. Here, ${ϵ}_b=8.9{ϵ}_0$.

Download Full Size | PDF

 

Fig. 9. Magnitude of $g_P^S(k,{\overline{k}}_i;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ as a function of ${\overline{k}}_i$ at $f_n=0.26$: (a) ${ϵ}_b=8.9{ϵ}_0$ and (b) ${ϵ}_b=8.9(1+0.11i){ϵ}_0$.

Download Full Size | PDF

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 ${ϵ}_b=8.9{ϵ}_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 ${ϵ}_b=8.9(1+0.11i){ϵ}_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.

 

Fig. 10. Spatial variations (magnitudes) of $g^S(k;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ in a lossless background of ${ϵ}_b=8.9{ϵ}_0$: (a) $f_n=0.1$, (b) $f_n=0.2$, and (c) $f_n=0.4$.

Download Full Size | PDF

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 ${ϵ}_b=8.9{ϵ}_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 ${ϵ}_b=8.9(1+0.11i){ϵ}_0$. This yields a complex $k$, and avoids the poles in the modal expansion. The agreement for the lossy case is better.

 

Fig. 11. $g^S(k;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ as a function of the normalized frequency: (a) ${ϵ}_b=8.9{ϵ}_0$ and (b) ${ϵ}_b=8.9(1+0.11i){ϵ}_0$.

Download Full Size | PDF

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.

 

Fig. 12. Spatial variations (magnitudes) of $g^S(k;\overline{\rho },{\overline{\rho }}^{\prime \prime })$ in a lossy background of ${ϵ}_b=8.9(1+0.11i){ϵ}_0$: (a) $f_n=0.1$, (b) $f_n=0.2$, and (c) $f_n=0.4$.

Download Full Size | PDF

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