In this work, we propose a numerical solver combining the spectral element - boundary integral (SEBI) method with the periodic layered medium dyadic Green’s function. The periodic layered medium dyadic Green’s function is formulated under matrix representation. The surface integral equations (SIEs) are then implemented as the radiation boundary condition to truncate the top and bottom computation domain. After describing the interior computation domain with the vector wave equations, and treating the lateral boundaries with Bloch periodic boundary conditions, the whole computation domains are discretized with mixed-order Gauss- Lobatto-Legendre basis functions in the SEBI method. This method avoids the discretization of the top and bottom layered media, so it can be much more efficient than conventional methods. Numerical results validate the proposed solver with fast convergence throughout the whole computation domain and good performance for typical multiscale nano-optical applications.
© 2017 Optical Society of America
The past decades witnessed the rapid developments in nano-optics, where various complex nano-scale designs and technologies evolve from laboratory prototypes to practical engineerings. Extreme ultra-violet (EUV) lithography, nano antennas and plasmonic metamaterials are the representatives of this large category of nano-scale technologies highlighted in the industrials [1,2]. Accompanying the engineering level applications, their design and fabrication increasingly address the importance of quantified analyses [3–5]. Various computational methods have been proposed for the analysis of nano-optical scattering problems. However, at optical spectrum, the target problem is usually electrically large, especially for 3-dimensional full-wave analysis. The versatile structure of the lithography pattern, optical connection and the background masks and substrate further constitute a multiscale problem. Conventional computational methods, such as the finite element method and the finite-difference time-domain method, have been under intensive study for their flexibility and robustness in this area . Unfortunately, these approaches have many difficulties when handling multiscale problems. Their slow convergence also requires a particular high sampling density inside the computation domain. Solving the final linear equation system usually consumes huge computational resources .
The spectral element - boundary integral (SEBI) method is among the most promising methods for electromagnetic scattering analyses [7,8]. Inherited from the philosophy of the finite element - boundary integral (FEBI) method, its computational framework incorporates the exact radiation boundary condition with only a minor increase of the knowns, while its strong flexibility keeps the SEBI a competitive solution for a wide category of problems [9–11]. One of the major drawbacks of the conventional FEBI is its slow convergence. In many cases, a minimum 20 points per wavelength (PPW) sampling density is required for an acceptable accuracy. As an alternative, the SEBI framework implements the mixed-order Gauss-Lobatto-Legendre (GLL) sampling strategy, so that fast convergence can be achieved throughout the computation domain. Intensive investigations have been reported for the application of SEBI in microwave and optics focusing on the electromagnetic scattering from objects embedding in homogeneous or non-periodic environment [11–13]. In spite of its efficiency and robustness, electrically small periodic nano structures fabricated in an electrically large planar layered medium, however, consists of the most common scenario in optical designs. In this work, we focus on the solutions for the periodic objects within layered medium background. For this category of problems, the conventional finite element boundary integral method requires discretization of the whole multi-scale model, resulting in a unnecessarily large number of unknowns. The surface integral equation (SIE) and volume integral equation (VIE) are also applicable approaches for the layered medium background problems. However, the SIE can only deal with homogeneous scatterers, while the VIE usually arrives at a typically large and dense linear equation system [14,15].
In this paper, we extend the SEBI framework to the study of periodic heterogeneous scatterers in planar layered medium background. Under the SEBI framework, the interior computation domain is formulated with the weak-form vector wave equation, while the radiation boundary is implemented with a set of surface integral equations. To describe the general periodic scatterers embedded in a layered medium, the SIE needs to be formulated with dyadic Green’s functions. In addition, the periodic layered medium Green’s function (PLMGF) must be carefully formulated to deliver a satisfying convergence rate together with a friendly implementation interface. With the combined SEBI-PLMGF method, the computation domain completely excludes the layered medium background, and an exponential convergence is expected for the solver because of the mixed-order GLL discretization for both the spectral element method (SEM) domain and the boundary integral domain.
The rest of the paper is organized as follows: Section 2 provides the overall framework of the SEBI for general periodic scatterers in a multilayered medium. Section 3 discusses the formulation and implementation of PLMGF within the context of SEBI. The treatment of primary field terms and secondary field terms are implemented separately according to their analytical characteristics. For the analysis of periodic structures, the SIE only serves as the radiation boundary condition to truncate the top and bottom of the computation domain. The lateral boundaries are further related with the Floquet theorem. A short note about combining the Bloch periodic boundary condition with SEBI is included in Section 4. Finally, several typical numerical results are provided in Section 5 to validate the proposed solver. The computational resources consumed by different approaches are also illustrated within each example.
2. Periodic boundary integral spectral element method for layered background
In this section, we first derive the overall framework of the SEBI. The overall computational framework is illustrated through Fig. 1.
Inherited from the conventional FEBI method, the SEBI-PLMGF method truncates the open boundary surface of the computation domain with a set of SIEs, while the interior computational domain is interpolated with a special category of local basis functions. In order to effectively model the multiscale nano-optical scattering problem, the SIEs are formulated with the periodic layered medium dyadic Green’s function, thus, the unknowns for the wave equation describing the background formation can be compressed into a set of equivalent surface integral equations. The region without scatterers therefore can be truncated from the final computation domain. Moreover, the Floquet theorem is further enforced as a set of periodic boundary conditions at the final linear equation system level, hence, the computation domain only cover one principal cell of the whole periodic structure.
Under the SEBI framework, the electric field in the interior computation domain is governed by the vector Helmholtz equation. Following the Galerkin’s method, the weak-form wave equation can be obtained as 
For the electromagnetic scattering problem with dual periodicity in the â1 and â2 directions, the radiation boundary consists of two separate surfaces. Without loss of generality, we assume the scatterer is periodic in x and y directions. Therefore, the radiation boundary consists of the top and bottom surfaces of the computation domain, both parallel to the XOY plane with no mutual coupling. For each of these surfaces, a set of electric field integral equation (EFIE) and magnetic field integral equation (MFIE) are combined with Eq. (1) to describe the wave propagation in the exterior domain:Fig. 1, each of the top and bottom boundary integral surfaces are governed by a set of EFIE and MFIE equations, separately.
There are generally four types of integral operators associated with the SIE: ℒE, 𝒦E, ℒH, 𝒦H,
Equation (2) (the EFIE) is directly applied with Galerkin’s method to arrive at its weak form:Eqs. (8) – (9) include both the top and the bottom radiation boundary surfaces.
Combining Eqs. (8) – (9) with proper discretization strategy, the electromagnetic scattering problem for periodic scatterers embedded in layered background can be solved without interior resonance concern [16,17]. The final matrix-form linear equation system will be derived in detail in Section 4 after the systematic discussion of the PLMGF and the periodic boundary conditions’ implementation.
3. Periodic layered medium Green’s function
The necessity of introducing the PLMGF originates from the SIE-based radiation boundary condition. The PLMGF here can be intuitively considered as the impulse response of a set of periodic electric/ magnetic point sources embedded in the layered background medium. For arbitrary electric and magnetic current sources, the governing equation of the field response can be conveniently expressed with the dyadic magnetic-type (ḠP,m) and electric-type (ḠP,e) periodic Green’s functions. In the remaining discussions, the subscript m and n denotes the layer index of the source point and observation point, respectively; p, q are the indices of the double summation; km is the wave number in the m-th layer.
The electric-type and magnetic-type PLMGF can be equivalently expressed with their spectral domain components or spatial domain components:Eqs. (4) and (7). The terms associated with ḠP,m and ∇ × ḠP,m can be easily obtained through the duality principle [14, 18]. For simplicity, the double summation operator is denoted as for the remaining of this article.
In general, for one-dimensional piecewise constant materials, the non-periodic layered medium Green’s function can be determined analytically in the spectral domain, and the spatial domain counterpart needs to be obtained through calculating the Sommerfeld integrals, which is usually a costly process. For periodic layered medium Green’s function, although closed formulation usually cannot be obtain in either domain, the Poisson’s summation theorem bridges the implementation roadmap in the spectral domain and the spatial domain. The double summation in spatial domain can be equivalently transformed into that in spectral domain.
It might be intuitive that computing the spectral domain summation is considerably superior to realizing the spatial domain formulation, where every single term in the series represents a shifted non-periodic dyadic layered medium Green’s function. Each term’s preparation requires a Fourier transform of its closed-form spectral domain counterpart. However, from the implementation point of view, the PLMGF is highly singular within the context of weak-form SIEs. Singularity issues are always present under pure spectral or spatial domain formulation. Conversely, the spectral domain and spatial domain representation exhibits distinct analytical properties. Switching between these two domains is essential to convert the generalized function into integrable express, and thus extract the singularities.
The same concern also holds for the primary terms and the secondary terms of the PLMGF for their distinct analytical properties. For this reason, we always split the Green’s function related SEBI terms into primary and secondary components, and formulate them separately. The singularity issues therefore can be extracted from the series summation and avoid massive unnecessary computation.
3.2. Matrix representation of the PLMGF’s primary terms
The primary term of the dyadic Green’s function is defined as the point source’s direct contribution in the source layer without any reflections from other layers. Therefore, the PLMGF’s primary terms are the same as the dyadic Green’s function for a homogeneous medium, and is non-zeros only when the observer and the source are within the same medium layer.
3.2.1. Primary term of ḠP,e
The primary term of ḠP,e in Eq. (4) is
The hyper-singularity in the above expression originates from the ∇∇ operator. Fortunately, within the context of SEBI, the spatial derivative can be shifted to the testing and basis functions. The primary components of the ḠP,e related terms can be express in matrix representation:
The integral kernel in Eq. (13) converges linearly with p, q under the spatial domain realization. As an alternative, applying the Poisson’s summation theorem, it can be equivalently expressed as
3.2.2. Primary term of ∇ × ḠP,e
The primary term of ∇ × ḠP,e in Eq. (7) is
The first attempt to accelerate the convergence of the above is transforming it into a spectral domain series. Through the Poisson’s summation theorem, Eq. (15) reaches at
Equation (16) exhibits exponential convergence when z ≠ z′. However, when z = z′ the above is a typical scenario where singularity issue arises: ∂ze−jkmz·|z−z′| results in a generalized function. As a result, the primary term of ∇ × ḠP,e needs to be formulated separately for two situations. When z ≠ z′, computing through Eq. (16) with exponential convergence; when z = z′, stay with Eq. (15) so that the singularity can be extracted. To guarantee the robustness of the overall solution, a numerical criteria should be given for switching between Eqs. (16) and (15). As one applicable approach, it can be set by the residual value of km · |z − z′| with respect to the specific precision adopted in implementation.
Although the above are valid in general, the surface integral kernel is still strongly singular (1/R2-type) when z = z′. However, recall that all the formulation of the PLMGF is prepared under the assumption that the scatterer is periodic in the x and y directions. Within the periodic SEBI framework, the scattering boundary surface is parallel to the XOY plane, i.e. n̂ = ±ẑ and z = z′. The xx, xy, yx, yy tensor components of ∇ × ḠP,e are zero at z = z′, thus considerably simplify the final implementation. The primary components of the ∇ × ḠP,e related terms can be expressed in matrix representation as
3.3. Matrix representation of the PLMGF’s secondary Terms
The secondary term of the PLMGF is defined to account for the transmission and reflection of the layered medium.
3.3.1. Secondary term of ḠP,e
The matrix representation of ḠP,e ’s secondary term is
The above spatial domain summation expressions of and directly originate from the definition of periodic Green’s function. They can be equivalently written as
Different from the PLMGF’s primary field terms, the summation kernel and its spatial derivatives of Eqs. (22) – (23) cannot be expressed in closed form. To avoid repeatedly calculating the inverse transformation or Sommerfeld integral, the above can be transformed into spectral domain realization:
3.3.2. Secondary term of ∇ × ḠP,e
For ∇ × ḠP,e in Eq. (7), matrix representation of the secondary term is
Similarly to the secondary term of the ḠP,e, all components inside the double summation operator are not in closed form. Fortunately, the above can be equivalently transformed into spectral domain realization:
It should be noticed that the secondary terms of the PLMGF converges exponentially if the source point and the observation point are not exactly on the same layered medium interface. This condition can be comfortably satisfied in most periodic SEBI applications. Even for the case z = z′ = zm, where zm denotes the z-coordinate of the m-th layered medium interface, the singularity still can be efficiently extracted under the matrix friendly formulation of the PLMGF .
4. Discretization and Bloch periodic boundary condition
Different from the philosophy of realizing the radiation boundary condition in the SEBI framework, the enforcement of the periodic boundary condition is derived under the spectral element method framework for boundary value problems. Before enforcing the Bloch boundary conditions to the lateral surfaces of the computation domain, we first discretize the unknowns in the weak-form wave equations and the SIE’s.
One of the major advantages of the SEBI or SEM originates from the GLL interpolation basis function. Compared with the quadratic convergence of the conventional finite element method , the SEM exhibits exponential convergence with respect to the increase of the order of interpolation polynomials. Discretizing both the SEM and MoM domain with the GLL basis function, Eqs. (8) – (9) arrive at the linear equation system:
Under the proposed SEBI framework for periodic scatterers embedded in layered medium, the top and bottom radiation boundaries have been realized by the surface integral equations with the PLMGF. The lateral boundaries of the computation domain are truncated with the Bloch periodic boundary condition:Eqs. (33) – (34) are only enforced on E’s tangential components on the periodic boundary surfaces of the whole computation domain, and H’s tangential components on the boundary edges of the top and bottom surfaces.
Equation (32) can be more compactly written as
To apply the Bloch boundary condition to the linear equation system, we first classify vector u according to each unknown’s location:
The Floquet theorem implies that
Equations (37) – (41) relate the corresponding unknown pairs on the PBC boundaries, so that the number of degree of freedom (DoF) can be reduced. The enforcement of Eqs. (37) – (41) consist of two steps: (a) remove the dependent testing functions in the weak-form equations (b) reduce the dimension of unknowns . This process can be elegantly formulated by introducing the Bloch phase matrix
The final linear equation system can be obtained as
All unknowns of the computation domain can be obtained by solving Eq. (43).
5. Numerical results
In this section, we show three typical numerical examples to validate the proposed solver and demonstrate its computational performance. First, the PLMGF is separately validated by comparing with the field results obtained from the finite element method. In the next example, we introduce a sub-wavelength nano scale example to validate the SEBI-PLMGF solver. Then, the fast convergence advantage of the SEBI-PLMGF is illustrated for three typical scenarios. Finally, the proposed solver is extended to electrically large multiscale extreme ultraviolet (EUV) lithography simulation, which is extremely challenging for conventional methods.
5.1. Case 1: Periodic dipole array in layered medium
In order to validate the PMLGF, the first numerical example is designed as an array of tilted electric dipoles placed in a three-layer medium. The material parameters of each layer is (from top to bottom): (1) ∊r = 1, μr = 1.5; (2) ∊r = 2, μr = 1; (3) ∊r = 3, μr = 1.5. The layered medium interfaces are at z = 0 nm and z = −30 nm. In the principal cell, an electric dipole is placed at (35, 35, −10) nm with the polarization of (θ = 30°, ϕ = 20°), and frequency of 7.5 × 1014 Hz. The x − and y−periodicities are 70 nm and 70 nm, respectively.
The electric and magnetic fields at the observation line (0, 50, −20) nm ∼ (70, 50, −20) nm are obtained through both the pure PLMGF package and the commercial software COMSOL, which utilizes the finite element method. The results of PLMGF are validated with the reference data in Fig. 2.
5.2. Case 2: Multi-layer nano structure
After validating the PLMGF, the second case is designed to validate the proposed SEBI-PLMGF solver. Figure 3 shows the geometric structure of the problem. The lengths of the principal cell in the x and y directions are 100 nm and 100 nm, respectively. The nano particle above the layered medium is of the size 50 nm ×50 nm ×50 nm, and the whole structure is under TMy with incident angle θ = 30°, ϕ = 20° and wavelength of 400 nm.
The SEBI-PLMGF solver and the conventional SEBI-homogeneous Green’s function solver are both used to simulate the same structure in Fig. 3. With the SEBI-PLMGF solver, the computation domain only contains the nano particle section. On the contrary, the SEBI-homogeneous Green’s function solver needs to include the whole structure inside the computation domain, which considerably increases the number of unknowns.
Observing on the plane z = 672.5 nm, the electric field distribution is shown in Fig. 4(a). Furthermore, the electric fields along the middle lines parallel to x axis and y axis are obtained from both solvers. From Fig. 4(b), a satisfying match between two sets of results can be observed. On an Intel 4-th generation i7 CPU, the SEBI-PLMGF approach takes 27s and 1.2GB memory to solve the 13,872 unknowns, while the SEBI-homogeneous Green’s function approach takes 59s and 4.5GB memory to solve the 72,800 unknowns. The processing work of meshing and discretization is not included in the time cost.
5.3. Case 3: p-refinement study
The major advantages of the proposed SEBI-PLMGF method includes two aspects. First, the layered medium background can be excluded from the computation domain. Only the fine structures of the principal cell need to be discretized. Second, the fast convergent GLL basis function is used to discretized both the interior computation domain and the boundary surfaces. In this section, several numerical examples are shown for the convergence performance of the SEBI-PLMGF.
As a commonly encountered issue in the conventional numerical methods, such as the finite element method, the finite-difference time-domain method, and the method of moments, the field singularities caused by the scatters’ sharp corner usually considerably affects the convergence of the numerical results. To address this concern, the refinement study is performed on three categories of problems, separately.
The scattering structures that impact the convergence performance of the SEBI-PLMGF solver can be considered as three categories in general: the computation domain has geometries with (a) no field singularities, (b) moderate singularities, and (c) significant field singularities. Figure 5 provides a sketch of the first two scenarios, while the previous multi-layer nano structure case, is a good representative for the third situation.
The p-refinement studies for all three situations are provided in Fig. 6. The A, B and C labels in the legends represent the computation model in Fig. 5(a), Fig. 5(b) and Fig. 3, respectively. As expected, the field singularities exhibit an impact on the convergence performance at various levels. Usually singular fields will cause the convergence to be slower than exponential. The numerical results, however, converge on an exponential scale in all cases.
It should also be noted that the numerical accuracy is not purely determined by the GLL basis functions’ order. Numerous factors, such as the mesh quality and the final linear equation system’s condition number, also have an non-negligible contribution for a well-behaved outcome.
5.4. Case 4: 3-dimensional extreme ultraviolet (EUV) lithography
With the above validation and performance discussion, we extend the proposed solver to a more realistic example. The EUV lithography is among the most typical demands for large scale numerical analysis of scatterers embedded in stratified structures. The overall structure is usually electrically large as each principal cells extents to dozens of wavelength in each direction. On the other hand, the actual size of the lithography pattern only occupies a small fraction of the whole structure. This category of large and multiscale problems are challenging for conventional numerical solvers to reach a satisfying accuracy with reasonable resources costs.
A EUV lithography model of typical principal cell size is designed for an illustration. Figure 7 provides a systematic sketch of the lithography model. Here we are etching a periodic pattern of DU embedded in multilayer structures. From top going downward, the materials and geometries are given as follows:
- Vacuum layer (infinitely thick)
- Top dielectric lens (5 bilayers, upper layer 3 nm-thick with ∊r = 0.998; lower layer 2 nm-thick with ∊r = 0.854)
- Lithography filling background (thickness: 10 nm, ∊r = 0.958)
- Absorber capping (thickness: 12 nm, ∊r = 0.7497 − j0.0296)
- Lithography pattern (thickness: 27 nm, ∊r = 0.856 − j0.0807)
- Multilayer capping (thickness: 2.5 nm, ∊r = 0.75 − j0.0296)
- Bragg reflector (40 silicon-molybdenum bilayers: Si layer thickness, 4.17 nm with ∊r = 0.998 − j0.00363; Mo layer thickness, 2.78 nm with ∊r = 0.854 − j0.0119)
- Substrate (silicon, infinitely thick)
Observing at one wavelength above the nano particle pattern, the scattered electric field distribution is shown in Fig. 8. The computation uses 15 min with 20 GB memory to solve the 65,790 unknowns for a converged result. Compared with Case 2, the increased computational resources demand is caused by the lower sparsity ratio of the final equation system. The relatively large-scale boundary integral surface results in a dense MoM block matrices with high dimension. The conventional FEBI method, however, is expected to require over 20 million unknowns for the same problem. Solving the massive equation system with the even larger MoM block matrices require orders of magnitude additional resources.
A mixed-order spectral element - boundary integral method for periodic scatterers embedded in a layered medium is systematically developed in this paper. Within the proposed SEBI-PLMGF framework, the minimum computation domain only contains the principal cell of the periodic scatterer. The layered medium background is described by a set of surface integral equations with the periodic layered medium dyadic Green’s functions, while the periodic boundaries are truncated with the Floquet theorem. The GLL polynomials based spectral edge basis functions are used to discretize the whole computation domain. The matrix representation of the dyadic PLMGF is also systematically discussed. The primary and secondary terms of the PLMGF exhibits distinct analytical properties. By formulating them separately, the strong singularity of the PLMGF can be considerably relieved. The theoretical foundation for combining the Floquet theorem with the SEM is also revisited in the paper. Several numerical examples are presented to validate the formulation. The proposed framework is expected to be an efficient solution for the electromagnetic analysis in a broad class of applications, such as lithography, solar cell, metamaterial, and optoelectronic designs.
References and links
1. B. Wu and A. Kumar, “Extreme ultraviolet lithography: A review,” J. Vac. Sci. Technol. 25, 1743–1761 (2007). [CrossRef]
3. J. Pomplun, S. Burger, L. Zschiedrich, and F. Schmidt, “Adaptive finite element method for simulation of optical nano structures,” physica status solidi (b) 244, 3419–3434 (2007). [CrossRef]
4. H. Lin, M. F. Pantoja, L. D. Angulo, J. Alvarez, R. G. Martin, and S. G. Garcia, “Fdtd modeling of graphene devices using complex conjugate dispersion material model,” IEEE Microw. Compon. Lett. 22, 612–614 (2012). [CrossRef]
5. I. Ahmed, E. H. Khoo, and E. Li, “Efficient modeling and simulation of graphene devices with the lod-fdtd method,” IEEE Microw. Compon. Lett. 23, 306–308 (2013). [CrossRef]
6. K. Saitoh and M. Koshiba, “Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. 19, 405 (2001). [CrossRef]
7. J. Niu, M. Luo, Y. Fang, and Q. H. Liu, “Boundary integral spectral element method analyses of extreme ultraviolet multilayer defects,” J. Opt. Soc. Am. A 31, 2203–2209 (2014). [CrossRef]
8. J.-H. Lee, T. Xiao, and Q. H. Liu, “A 3D spectral-element method using mixed-order curl conforming vector basis functions for electromagnetic fields,” IEEE Trans. Microw. Theory Techn. 54, 437–444 (2006). [CrossRef]
9. M. Luo, Y. Lin, and Q. H. Liu, “Spectral methods and domain decomposition for nanophotonic applications,” Proc. IEEE 101, 473–483 (2013). [CrossRef]
10. J. Niu, M. Luo, J. Zhu, and Q. H. Liu, “Enhanced plasmonic light absorption engineering of graphene: simulation by boundary-integral spectral element method,” Opt. Express 23, 4539–4551 (2015). [CrossRef] [PubMed]
11. J. Niu, M. Luo, and Q. H. Liu, “Enhancement of graphene’s third-harmonic generation with localized surface plasmon resonance under optical/electro-optic Kerr effects,” J. Opt. Soc. Am. B 33, 615–621 (2016). [CrossRef]
12. J. H. Lee, J. Chen, and Q. H. Liu, “A 3D discontinuous spectral element time-domain method for Maxwell’s equations,” IEEE Trans. Antennas Propag. 57, 2666–2674 (2009). [CrossRef]
13. Y. Ren, W. F. Huang, J. Niu, and Q. H. Liu, “A hybrid solver based on domain decomposition method for the composite scattering in layered medium,” IEEE Antennas Wireless Propag. Lett. 16, 420–423 (2017). [CrossRef]
14. Y. P. Chen, W. C. Chew, and L. Jiang, “A new Green’s function formulation for modeling homogeneous objects in layered medium,” IEEE Trans. Antennas Propag. 60, 4766–4776 (2012). [CrossRef]
15. Y. Ren, Q. H. Liu, and Y. P. Chen, “A hybrid fem/mom method for 3-d electromagnetic scattering in layered medium,” IEEE Trans. Antennas Propag. 64, 3487–3495 (2016). [CrossRef]
16. B. Shanker, A. A. Ergin, K. Aygun, and E. Michielssen, “Analysis of transient electromagnetic scattering from closed surfaces using a combined field integral equation,” IEEE Trans. Antennas Propag. 48, 1064–1074 (2000). [CrossRef]
17. J. D. Collins, J. M. Jin, and J. L. Volakis, “Eliminating interior resonances in finite element-boundary integral methods for scattering,” IEEE Trans. Antennas Propag. 40, 1583–1585 (1992). [CrossRef]
18. W. C. Chew, J. L. Xiong, and M. A. Saville, “A matrix-friendly formulation of layered medium Green’s function,” IIEEE Trans. Antennas Propag. Lett. 5, 490–494 (2006). [CrossRef]
19. K. Chen, J. Song, and T. Kamgaing, “Accurate and efficient computation of layered medium doubly periodic Green’s function in matrix-friendly formulation,” IEEE Trans. Antennas Propag. 63, 809–813 (2015). [CrossRef]
20. J. Jin, The Finite Element Method in Electromagnetics (Wiley-IEEE Press, 2014), 3rd ed.
21. A. Nicolet, S. Guenneau, C. Geuzaine, and F. Zolla, “Modelling of electromagnetic waves in periodic media with finite elements,” Journal of Computational and Applied Mathematics 168, 321–329 (2004). Selected Papers from the Second International Conference on Advanced Computational Methods in Engineering (ACOMEN 2002). [CrossRef]