1. Introduction

The band solutions of waves in periodic structures are important in physics and engineering and in the design of photonic, electronic, acoustic, microwave and millimeter wave devices such as that in photonic crystals and metamaterials. The common approaches to calculate the bands include the planewave method [1–7], and the Korringa Kohn Rostoker (KKR) method [8, 9]. The planewave 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 planewave method, the permittivity (or potential) profile is expanded in Fourier series. For the case of abrupt profiles of large contrasts between background and scatterer, many terms are required in the Fourier series expansion. An artificially smeared dielectric function was also employed [5]. The convergence issue of the planewave method has been reported [5, 7]. The KKR method is a multiple scattering approach [10, 11] by forming the surface integral equations with periodic Green’s function and the equation is solved by cylindrical or spherical waves expansions. The method is only convenient for scatterers with shapes of separable geometries such as circular cylinders or spheres. In the KKR/multiple scattering method, the eigenvalue problem is nonlinear requiring the root seeking procedure. One needs to repeat the non-linear root seeking one by one to find multiple band solutions. Numerical approaches such as the finite difference time domain method (FDTD) [12–13] and finite element method (FEM) [14–16] have also been used.

In this paper, the periodic Green’s function is used to formulate the dual surface integral equations. The periodic Green’s function has slow convergence. If the periodic Green’s function is used for the surface integral equation to calculate the eigenvalues, the equation is nonlinear and an iterative search needs to be performed with one band at a time. Recently, the Broadband Green’s Function with Low wavenumber extraction (BBGFL) [17–20] has been applied to wave propagation in waveguide/cavity of arbitrary shape. By using a single low wavenumber extraction, the convergence of the modal expansion of the Green’s function is accelerated. Also the singularity of the Green’s function has been extracted. The method has been shown to be efficient for broadband simulations of wave propagation in waveguides/cavity. It was noted that the expression of modal expansion of Green’s function in a waveguide is similar to the Floquet expansion of Green’s function in a periodic structure [21]. In a recent paper [21] we adapted the BBGFL to calculate band diagrams in periodic structures where the scatterers obey Dirichlet boundary conditions. We use a low wavenumber extraction to accelerate the convergence of the periodic Green’s function. The BBGFL is used to formulate the surface integral equations. Next applying the Method of Moments (MoM) gives a linear eigenvalue equation that gives all the multi-band solutions simultaneously for a point in the first Brillouin zone. We label this as “broadband simulations” as the multi-band solutions are calculated simultaneously rather than searching the band solution one at a time. The BBGFL method was shown to be accurate and computationally efficient. The choice of the low wavenumber was shown to be robust as the choice of the low wavenumber can be quite arbitrary. The application of MoM makes the method applicable to arbitrary shapes. It is noted that the Boundary Integral Resonant Mode Expansion (BIRME) method [22] was used to find the modes of an arbitrary shaped waveguide with PEC or perfect magnetic conductor (PMC) boundary conditions, using the Greens functions of a rectangular waveguide or circular waveguide. In BIRME, a DC extraction is used for the Greens function. In this paper, we calculate band diagrams of periodic structures with general shapes of the unit cells and with dielectric scatterers embedded. Also, a low wavenumber extraction is used and the choice of the low wavenumbers is shown to be robust.

This paper is an extension of the recent paper of BBGFL on perfect electric conductor [21] to periodic dielectric scatterers. The case of dielectric scatterers have much wider applications in devices of photonic crystals and metamaterials. For the case of dielectric scatterers, the dual surface integral equations are derived. A distinct feature in the formulation is that we use the periodic Green’s function of the scatterers, in addition to using the periodic Green’s function of the background. The application of the broadband periodic Green’s function separates out the wavenumber dependence of Green’s function, and allows the conversion of the non-linear root searching problem into a linear eigenvalue analysis problem. The low wavenumber extractions are applied to both periodic Green’s functions for fast Floquet mode convergence and extraction of Green’s function singularities. The use of MoM makes the method applicable to scatterers of arbitrary shape. In section 2, we formulate the dual surface integral equation using the extinction theorems and the periodic Green’s functions. In Section 3, we apply BBGFL to the surface integral equations and derive the linear eigenvalue problem. In section 4, numerical results are illustrated for the band solution of circular air voids forming a honeycomb structure. The results are in good agreement with the plane wave method and the KKR method. We also describe the quick rejection of spurious modes by using extinction theorems. In Section 5, we use the lowest band around the Γ point to calculate the low frequency dispersion curves, the effective permittivity, and the effective propagation constants, all of which have wide applications in devices of metamaterials.

2. Extinction theorem and surface integral equations

Consider a periodic array in the xy plane of 2D dielectric scatterers embedded in a 2D lattice as illustrated in Figure 1. Let the primitive lattice vectors be ${\overline{a}}_1$ and ${\overline{a}}_2$ with lattice constant a and the primitive cell area ${\mathrm{\Omega }}_0=\left|{\overline{a}}_1\times {\overline{a}}_2\right|$. The scatterer has dielectric constant ε_rp and the background host ε_rb. Let S_pq denote the boundary of the scatterer in the (p,q)-th cell and A_pq is the domain of the scatterer in the (p,q) cell bounded by S_pq.

 

Fig. 1 Geometry of the 2D scattering problem in 2D lattice. 2D dielectric scatterers form a periodic array with primitive lattice vectors ā₁ and ā₂, and lattice constant a. The primitive cell area Ω₀ = |ā₁ × ā₂|. The (0, 0)-th scatterer has boundary S₀₀ and cross section A₀₀. The scatterer and background has dielectric constant ε_rp and ε_rb, respectively.

Download Full Size | PDF

We denote the field outside the scatterer as ψ and inside as ψ₁, where ψ represents E_z for TMz polarization and H_z for TEz polarization. TMz means that the magnetic field is perpendicular to the z direction while TEz means that the electric field is perpendicular to the z direction. The extinction theorems state that

In (1a) and (1b) $g(\overline{\rho },{\overline{\rho }}^{\prime })=\frac{i}{4}{H}_0^{(1)}(k\rho )$ and $g_1(\overline{\rho },{\overline{\rho }}^{\prime })=\frac{i}{4}{H}_0^{(1)}(k_1\rho )$ are the 2D free space Green’s functions in the the background and in the scatterer, respectively. The wavenumbers are $k=\frac{\omega }{c}\sqrt{{\epsilon }_{rb}}$, and $k_1=\frac{\omega }{c}\sqrt{{\epsilon }_{rp}}$ in the background and in the scatterer, respectively, where ω is the angular frequency and c is the speed of light.

Applying the Bloch wave condition to Eq. (1a), we have the surface integral to be only over S₀₀ of the (0,0) cell. if $\overline{\rho }$ is inside any scatterer, then

Let ${\overline{k}}_i$ be a wavevector in the first irreducible Brillouin zone ${\overline{k}}_i={\beta }_1{\overline{b}}_1+{\beta }_2{\overline{b}}_2$, 0 ≤ β₁, Β₂ ≤ 1/2. ${\overline{b}}_1$ and ${\overline{b}}_2$ are the reciprocal lattice vectors. ${\overline{b}}_{{}_1}=2\pi \frac{{\overline{a}}_2\times \widehat{z}}{{\mathrm{\Omega }}_0}$, ${\overline{b}}_2=2\pi \frac{\widehat{z}\times {\overline{a}}_2}{{\mathrm{\Omega }}_0}$. Then in the spectral domain the periodic Green’s functions for, respectively, background and domain scatterer are

The surface unknowns are $\psi ({\overline{\rho }}^{\prime })$ and ${\psi }_1({\overline{\rho }}^{\prime })$, ${\widehat{n}}^{\prime }\cdot {\nabla }^{\prime }\psi ({\overline{\rho }}^{\prime })$ and ${\widehat{n}}^{\prime }\cdot {\nabla }^{\prime }{\psi }_1({\overline{\rho }}^{\prime })$. There are boundary conditions for the TM case and the TE case giving a net of two surface unknowns $\psi ({\overline{\rho }}^{\prime })$ and ${\widehat{n}}^{\prime }\cdot {\nabla }^{\prime }\psi ({\overline{\rho }}^{\prime })$ [10,23]. Let $\overline{\rho }$ approach the boundary of S₀₀, the coupled surface integral equations in $\psi ({\overline{\rho }}^{\prime })$ and ${\widehat{n}}^{\prime }\cdot {\nabla }^{\prime }\psi ({\overline{\rho }}^{\prime })$ are given as follows,

Applying MoM with pulse basis functions and point matching, we discretize S₀₀ into N patches. In matrix form, we obtain

3. BBGFL in surface integral equations

We want to find the wavenumbers that satisfy both the surface integral equations in Eq. (4). These wavenumbers are the band solutions, and the corresponding surface field ψ and $\widehat{n}\cdot \nabla \psi$ are the modal currents. This is an eigenvalue problem. The technique of BBGFL is described below.

3.1. Low wavenumber extraction in the Broadband Periodic Green’s function

In the BBGFL method, we choose a single low wavenumber k_L, which is quite arbitrary, to decompose $g_P(k,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$ into $g_P(k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$ and the difference $g_B(k,k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$. They shall be labelled as the low wavenumber part and the broadband part, respectively.

Using the spectral domain of $g_P(k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$ from (3a), we then have

Note that $g_P(k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$ is independent of wavenumber k, and the series in $g_B(k,k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$ converges as ${\left|k_{i\alpha }\right|}^{-4}$ instead of ${\left|k_{i\alpha }\right|}^{-2}$.

For the case of dielectric scatterer, we also decompose the periodic Green’s function of the scatterer $g_1{}_P(k,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$ into a low wavenumber part $g_1{}_P(k_{1L},{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$ with k_1L and a broadband part $g_1{}_B(k,k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$

We next express g_B in the following form

It follows that

Both g_B and ${\widehat{n}}^{\prime }\cdot {{\nabla }^{\prime }}_{gB}$ are smooth for arbitrary $\overline{\rho }-{\overline{\rho }}^{\prime }$ and they converge everywhere including the self-point of $\overline{\rho }={\overline{\rho }}^{\prime }$. Because of fast convergence, we will truncate the series with M Floquet modes. Note the two dimensional index in α. (In our numerical implementation, each index is from −10 to 10, so that M = (21)² = 441.)

With the decomposition $g_P(k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })+g_B(k,k_L,{\overline{k}}_i,\overline{\rho },{\overline{\rho }}^{\prime })$, the matrices $\overline{\overline{S}}$ and $\overline{\overline{L}}$ (Eq. 6) are also in two parts which are the low wavenumber part and the broadband part

$\overline{\overline{R}}$ and $\overline{\overline{Q}}$ are N × M matrices. The superscript ^† means Hermitian adjoint. The matrices $\overline{\overline{S}}(k_L),\overline{\overline{L}}(k_L),\overline{\overline{R}}(k_L)$ and $\overline{\overline{Q}}(k_L)$ do not have wavenumber dependence. The only wavenumber dependence lies in the factor $\lambda (k,k_L)$ of the M × M diagonal matrix $\overline{\overline{W}}(k_L,k)$

With the low wavenumber separation $\overline{\overline{S}}$ and $\overline{\overline{L}}$, the matrix form of the surface integral equation Eq. (5a) becomes

Using the surface integral equations with periodic Green’s functions of the scatterer $g_{1P}\left(k_1,{\overline{k}}_i;\overline{\rho },\overline{{\rho }^{\prime }}\right)$ which is also decomposed into a low wavenumber part and a broadband part, we have from Eq. (5b)

We also choose

In Appendix A, we describe efficient methods for calculating the matrix elements of $\overline{\overline{S}}$, $\overline{\overline{L}}$, ${\overline{\overline{S}}}^{(1)}$ and ${\overline{\overline{L}}}^{(1)}$, at low wavenumbers of k_L and k_1L.

3.2. Eigenvalue problem

We next convert the matrix equation Eq. (14) into an eigenvalue problem. Let

Then from Eq. (13a)

Thus

Similarly, with Eq. (16c–d, 15c)

Next, using the definitions of $\overline{b}$, $\overline{c}$, ${\overline{b}}^{(1)}$, ${\overline{c}}^{(1)}$, the matrix equation Eq. (14) becomes,

Next we express $\overline{p}$, $\overline{q}$ in $\overline{b}$, $\overline{c}$, ${\overline{b}}^{(1)}$, ${\overline{c}}^{(1)}$,

The matrix dimension of the eigenvalue problem in (23a) is 4M × 4M. In the above equation, only λ depends on wavenumber while all the other matrices are independent of wavenumber. Thus the eigenvalue problem is linear with all the eigenvalues and eigenvectors solved simultaneously giving the multi-band solutions.

Knowing the eigenvalues λ, the mode wavenumbers k² are obtained from the relation of $\lambda =\frac{1}{k^2-k_L^2}$, (Eq. (9c)). The modal surface currents distributions $\overline{p}$ and $\overline{q}$ are calculated from the eigenvectors through Eq. (20). Knowing $\overline{p}$ and $\overline{q}$, we can also compute the modal field distribution ψ and ψ₁ everywhere in the lattice by Eq. (2). The authenticity of the eigenmodes can also be checked by the extinction theorem (Eq. (2)) away from the boundaries.

4. Numerical results

We consider the triangular lattices with lattice constant a as shown in Figure 1. The dielectric background has permittivity ε_b, and the scatterers are circular air voids of radius b drilled in the background with ε_p = ε₀, where ε₀ is the free space permittivity. We consider two cases with case 1 b = 0.2a, ε_b = 8.9ε₀ and filling ratio of 14.5% and case 2 of b = 0.48a, ε_b = 12.25ε₀ and filling ratio of 83.6%. We calculate the solutions using BBGFL method. We also calculate the solutions using the plane wave method [4] and the KKR method [10] and compare the results of BBGFL with these two methods.

In Appendix A, we describe the calculations of g_P and g_1P at the respective single low wavenumbers of k_L and k_1L that include the calculations of D_n [10]. The lattice vectors are given by

The Γ, M, and K points are

We plot the band solutions with ${\overline{k}}_i={\beta }_1{\overline{b}}_1+{\beta }_2{\overline{b}}_2$ following Γ → M→ K → Γ.

The normalized frequency f_N (k) is defined as

In the following computations, we choose the low wavenumber k_L such that f_N (k_L) = 0.2. The choice of k_L is quite arbitrary as k_L values in a range will work. This makes BBGFL method robust. In the examples, the maximum number of two-dimensional (2D) Bloch waves in BBGFL is 441 corresponding to −10 ≤ m, n ≤ 10, so that the highest m, n index is 10. This leads to an eigenvalue problem of maximum dimension 4M = 1764.

Case 1: b = 0.2a, ε_b = 8.9ε₀, πb²/Ω₀ = 14.5%

The band diagram is plotted in Figure 2, with MoM discretization N = 80, and D_n (as defined in Eq. (27)) at order 4 to compute g_P and ${\widehat{n}}^{\prime }$ ∇′g_P at the single low wavenumber k_L. We compute solutions using the plane wave method with 1681 Bloch waves. Note that we only use 441 Bloch waves in the BBGFL method. The agreements between the two methods are excellent. In Figure 3, the surface modal currents, unnormalized, are plotted for the first few modes near Γ point with ${\overline{k}}_i=0.05{\overline{b}}_1$.

 

Fig. 2 Band diagram of the hexagonal structure with background dielectric constant of 8.9 and air voids of radius b = 0.2a. The results of BBGFL are shown by the solid curves for TMz polarization and by the dashed curves for the TEz polarization. The circles and crosses are the results of planewave method for the TMz and TEz polarizations, respectively.

Download Full Size | PDF

 

Fig. 3 Modal surface currents distribution near Γ point at ${\overline{k}}_i=0.05{\overline{b}}_1$ corresponding to the first few modes of the hexagonal structure with background dielectric constant of 8.9 and air voids of radius b = 0.2a. (a) TMz, surface electric currents $J_z\propto \partial \psi /\partial n$; (b) TMz, surface magnetic current M_t ∝ ψ; (c) TEz, surface magnetic currents $M_z\propto \partial \psi /\partial n$; (d) TEz, surface electric current J_t ∝ ψ. The corresponding normalized mode frequencies are 1: 0.0208, 2: 0.3736, 3: 0.3864 for TMz wave, and 1: 0.02224, 2: 0.3847, 3: 0.404 for TEz wave, respectively.

Download Full Size | PDF

In Table 1, the convergence of the lowest mode with respect to the number of Bloch waves used in BBGFL using different low wavenumber k_L are tabulated for TMz polarization with ${\overline{k}}_i=0.05{\overline{b}}_1$. It is noted that the choice of k_L is robust, and fewer Floquet modes are needed as one choses a lower k_L.

Table 1. The convergence of the lowest mode with respect to the number of Bloch waves used in BBGFL using different low wavenumber k_L. The results are tabulated for TMz polarization with ${\overline{k}}_i=0.05{\overline{b}}_1$, where 0.020798 is the first band solution.

View Table

Case 2: b = 0.48a, ε_b = 12.25ε₀, πb²/Ω₀ = 83.6%

The parameters for this case are the same as in [10]. The area ratio of scatterers to background is high at 83.6%. The dielectric constant ratio of ε_b/ε_p is as high as 12.25. In this case, the field varies more rapidly along the perimeter of the scatterer than the previous case and we choose the MoM discretization to be N = 180, and D_n (as defined in Eq. (27)) is selected at order 8 to compute g_P and $\widehat{n^{\prime }}\cdot {\nabla }^{\prime }{g}_P$ at the single low wavenumber k_L. The band diagram is plotted in Figure 4. For the planewave method we used 6561 Bloch waves which is significantly larger than M = 441 Floquet modes for BBGFL. The KKR solution of the first few modes at the points of M, and K are also given (Ref. [10] has the complete KKR results). The band gap between the third and fourth bands is noted. The agreements between the three methods are in general good. Note that the planewave method predicts slightly larger higher order modes than the BBGFL. The differences decrease as more Bloch waves are included in the planewave expansion. Note that the planewave method results shown are not exactly the same as the results reported in [10], possibly due to the different treatment of the Fourier series expansion of the dielectric function [4, 5]. In our implementation, we do not apply any smearing function to the permittivity profile ε(r), and directly compute the Fourier transform of ε⁻¹ ( $\overline{r}$). For the KKR method, more terms in the D_n series are needed as frequency increases. In Figure 5, the unnormalized surface modal currents are plotted for the first few modes near Γ point with ${\overline{k}}_i=0.05{\overline{b}}_1$. The modal currents tend to have more variations as the modal wavenumbers increase.

 

Fig. 4 Band diagram of the hexagonal structure with background dielectric constant of 12.25 and air voids of radius b = 0.48a. The results of BBGFL are shown by the solid curves for TMz polarization and by the dashed curves for the TEz polarization. The circles and crosses are the results of planewave method for the TMz and TEz polarizations, respectively. The triangles and squares are the results of the KKR method for the TMz and TEz polarizations, respectively.

Download Full Size | PDF

 

Fig. 5 Modal surface currents distribution near Γ point at ${\overline{k}}_i=0.05{\overline{b}}_1$ corresponding to the first few modes of the hexagonal structure with background dielectric constant of 12.25 and air voids of radius b = 0.48a. (a) TMz, surface electric currents $J_z\propto \partial \psi /\partial n$; (b) TMz, surface magnetic current M_t ∝ ψ; (c) TEz, surface magnetic currents $M_z\propto \partial \psi /\partial n$; (d) TEz, surface electric current J_t ∝ ψ. The corresponding normalized mode frequencies are 1: 0.03437, 2: 0.4425, 3: 0.6154 for TMz wave, and 1: 0.04145, 2: 0.7669, 3: 0.7710 for TEz wave, respectively.

Download Full Size | PDF

The eigenvalues of the BBGFL approach include spurious modes. All the spurious modes are quickly rejected based on the following simple calculations. The spurious modes do not satisfy the extinction theorem as given in Eq. (2). If we evaluate the field inside the scatterers using g_P, or evaluate the field outside the scatterers using g_1P, we get non-zero values for the spurious modes. These modes can be identified by several points calculations inside or outside the scattering using the surface integrals. Some of the spurious modes do not satisfy one of the two surface integral equations Eq. (4), and for these modes, their modal currents ψ and $\partial \psi /\partial n$ on the boundary are essentially zeros and trivial. For the rest of the spurious modes, they satisfy the surface integral equations Eq. (4), but do not obey the extinction theorem Eq. (2), and these modes are shown to be nearly invariant with respect to ${\overline{k}}_i$, with values close to the roots of J_n (kb), where J_n is the n-th order Bessel function.

The CPU time is recorded when running the BBGFL code in Mat lab R2014b to compute the band diagram on an HP ProDesk 600 G1 desktop with Intel Core i7-4790 CPU @ 3.60 GHz and 32 GB RAM. Since the computation for different ${\overline{k}}_i$ are independent, the CPU time for one ${\overline{k}}_i$ will be described. For case 2 with the larger N and higher order of D_n, it takes a total of ^~96 sec to complete the mode analysis at one ${\overline{k}}_i$, of which ^~11s is spent on the eigenvalue analysis. More than 80% CPU time is spent on computing the matrices $\overline{\overline{S}},\overline{\overline{L}},{\overline{\overline{S}}}^{(1)}$, and ${\overline{\overline{L}}}^{(1)}$ at the single low wavenumbers k_L and k_1L respectively for background periodic Green’s function and scatterer periodic Green’s function. The calculation of D (and $D_n^{(1)}$ by changing k to k₁) as defined in Appendix A to facilitate the computation of g_P and g_1P takes ^~24s, and the calculation of the low wavenumber matrices using D_n takes ^~27s on the boundary. Another ^~28s are needed for solutions away from the boundary as needed in checking the extinction theorem. The time recorded is for one polarization, and for the second polarization the added time is very small since only the eigenvalue problem is to be resolved. The approach is much more efficient than the KKR approach which is based on the evaluation of D_n at every frequency that needed to be used to search the band solution k, with one band at a time. The BBGFL can also be more efficient than the planewave method when the planewave method requires significantly larger number of Bloch waves as in the second case with large dielectric contrast and large filling ratio. Note that for the case of infinite contrast of perfect electric conductor (PEC), our previous paper [21] shows that BBGFL also works.

5. Low frequency dispersion relations, effective permittivity and propagation constants

In this section, we illustrate the low frequency dispersion relations which are useful for the design of devices in photonics and metamaterials.

In metamaterials, the scatterers and the lattice spacings are subwavelengths. Effective permittivities and effective propagation constants have been calculated for random media of dielectric mixtures [23–28]. In random medium, the positions of scatterers are random creating random phase in propagating waves. Every realization has different positions of scatterers although each realization has the same statistics and the field solutions of different realizations are different. The waves are decomposed into coherent waves and incoherent waves. The coherent wave has definite amplitude and phase while the incoherent waves have random amplitudes and phases giving speckle. In low frequency, the coherent waves dominate while at higher frequencies, the incoherent waves dominate. In random medium, the effective permittivity and the effective propagation constant are that of the coherent waves. On the other hand, the geometry in a periodic structure is deterministic with a single solution. The quasistatic method has been used to calculate the effective permittivity of a periodic medium [28]. It is interesting to note that at very low frequency, the effective permittivities as derived for random medium and for periodic medium are the same. In the following, we associate effective permittivity and effective propagation constant with the low frequency dispersion relation of the lowest band in the vicinity of the Γ point.

In Figures 6 and 7, we plot the ω − k lowest band near the Γ point for case a and case b, respectively. These will be labeled as low frequency dispersion curves. The dispersion relation curves are plotted for ${\overline{k}}_i$ moving along the line of ΓM, and along the line of ΓK in the first Brillouin zone, respectively. The points M and K represents the largest anisotropy in the first irreducible Brillouin zone. In Appendix B we derived the effective permittivity and effective propagation constants [23]. The corresponding dispersion curves are straight lines and are also plotted in Figures 6 and 7.

 

Fig. 6 Dispersion relationship of the hexagonal structure with background dielectric constant of 8.9 and air voids of radius b = 0.2a.

Download Full Size | PDF

 

Fig. 7 Dispersion relationship of the hexagonal structure with background dielectric constant of 12.25 and air voids of radius b = 0.48a.

Download Full Size | PDF

The low frequency dispersion relations of the lowest band are close to the effective permittivity curves at low frequency and deviate as $\left|{\overline{k}}_i\right|$ increases. At higher frequencies, towards the M and K points, they show significant departures when ${\overline{k}}_i$ are close to M and K point. Thus as frequency increases, the effective propagation constants and effective permittivity are dispersive and anisotropic. The dispersion relations are plotted for both TM polarization and TE polarization. For the same ω, TM wave in general has a slightly larger k. The difference between the two polarization increases as the permittivity contrast and the scatterer filling ratio increases.

The BBGFL method also works for infinite permittivity contrast as in the case of PEC [21]. In Figure 8, we plot the ω − k dispersion relation for the PEC scatterer case, with the PEC cylinder radius b = 0.2a, and the background permittivity ε_b = 8.9ε₀. The behavior of the TE wave is similar to that of the dielectric case. But for the TM wave, as indicated by the quasistatic mixing formula that the mixture does not have a finite quasistatic effective permittivity. The ω − k dispersion relation from the lowest band has similar behavior to that of the plasma [29]. The periodic structure behaves like a plasma for the TM wave that the wave could only propagate when its frequency is larger than the corresponding plasma frequency.

 

Fig. 8 Dispersion relationship of the hexagonal structure with background dielectric constant of 8.9 and PEC cylinders of radius b = 0.2a.

Download Full Size | PDF

6. Conclusions

In this paper we have extended the BBGFL approach, previously used for Dirichlet boundary conditions [21], to calculate band solutions and the low frequency dispersion relations of dielectric periodic structures in 2D problem with 2D periodicity. In the BBGFL approach, we solve the dual surface integral equations with MoM using the broadband periodic Green’s function of both the background and the scatterer. The Broadband periodic Green’s functions are fast convergent because of low wavenumber extraction. Because MoM is applied, the method is applicable to arbitrary scatterer shapes and arbitrary filling ratio and permittivity contrasts. The BBGFL approach has the form of a linear eigenvalue problem. The method is shown to be efficient and accurate. The choice of the low frequency wavenumber has been shown to be robust because MoM is applied. The method, in principle, is applicable to arbitrary scaterer shapes and arbitrary filling ratio and permittivity contrasts. We are studying these cases. We are also extending the method to 3D problems with 3D periodicity.