Abstract
A new, rigorous, field-based, seminumerical analysis method is presented to obtain the reflection and transmission coefficients of 2D planar periodic structures with arbitrarily shaped metallization patterns for both normal and oblique incidence conditions. It is useful for the analysis, design, and optimization of many single-layer and multilayer planar structures, such as frequency-selective surfaces (FSSs), artificial magnetic conductor (AMC) surfaces, electromagnetic bandgap (EBG) structures, some metamaterials and high-impedance surfaces. In this coupled-field expansion method (CFEM), the x- and y-components of the vector magnetic potential in each homogeneous region in a unit cell are expanded in terms of Bloch-Floquet modes and the solution to the coupled-field problem is formulated. The unique, analytical formulation presented here leads to a linear system with reasonably simple matrix elements. By cascading the matrices representing each interface, multilayer periodic structures are analyzed in a very flexible way. Being field based, CFEM does not require substrate Green's functions to analyze surfaces printed on dielectric substrates. The method was validated by analyzing one single-layer periodic surface (a printed AMC on a dielectric substrate) and one multilayer periodic surface (a circular polarizer) and comparing CFEM results with HFSS results.
1. Introduction
Periodic two-dimensional arrays of thin conducting patterns have been widely investigated as frequency-selective surfaces (FSS), artificial magnetic conductors (AMC), electromagnetic bandgap (EBG) structures, planar metamaterials, circular polarizers, and high-impedance surfaces. Fully numerical methods, such as the finite-element method (FEM) [1] and the finite-difference time-domain (FDTD) [2] method compute the field distribution in a unit cell completely numerically. The direct output of these methods, that is, numerical values of the fields at sampled points, does not indicate interesting information, such as which modes are propagating and which modes are evanescent, unless further processing is carried out.
The Bloch-Floquet theorem, which describes the periodic electromagnetic wave distribution in periodic structures, is an important starting point for many semianalytical methods. Periodic Green’s functions based on TE and TM polarized Bloch-Floquet expansions have been proposed and have been applied with the Method of Moments (MoM) to obtain the full-wave solution for 2D periodic arrays [3–8]. The electric or magnetic current distribution is the direct output of these methods. A single-layer FSS in free space was analyzed in [3] by combining TE and TM polarized Bloch-Floquet modes with the MoM. In [4], the MoM was extended for an FSS printed on a dielectric substrate. A set of entire-domain basis functions numerically calculated using the boundary integral resonant mode expansion (BI-RME) method was proposed for multilayer FSSs in [5]. In [6, 7], a Green’s function that takes into account the singularity at metal patch edges was exploited with entire-domain basis functions, and the singularity was treated analytically. A cascaded generalized scattering matrix was presented in [8] to efficiently solve multilayer problems. All these semianalytical methods solve an integral equation and compute the current distribution first. For substrate-based and printed structures, they require Green’s functions.
In this paper, we present a new, rigorous, field-based, seminumerical, full-wave analysis method for analyzing single-layer/multilayer planar periodic structures by analytically expanding the fields in each homogeneous region between interfaces. In the development of this method, the concepts were taken from established methods such as the MoM, Mode Matching Method, and Spectral Domain method. In this coupled-field expansion method (CFEM), the problem space is divided into two/several homogeneous regions which are coupled. The coupling occurs through aperture areas (e.g., apertures) in the planar interface between regions. The and components of the vector magnetic potential in each region is expanded using Bloch-Floquet modes. The unknown coefficients in the expansion are determined by analytically enforcing interface conditions over the entire area they are applicable. By characterizing each interface using a matrix and cascading the matrices, multilayer periodic structures are analyzed in a very flexible way. As we solve for fields instead of currents, substrate Green’s functions are not required to analyze multilayer periodic structures on dielectric substrates.
The theory for single-layer conducting patterns is described in Section 2.2. Then, CFEM is extended for multilayer case in Section 2.3. To verify the accuracy of this method, we present numerical results for a printed AMC on a dielectric substrate and a multilayer circular polarizer, and compare CFEM results with HFSS results in Section 3. The CFEM results agree very well with HFSS simulation results for normal and oblique incidence. The surface designed by this method has been successfully employed in the small EBG-resonator antennas.
2. Thoery
2.1. Formulation of the Problem
Due to translational periodicity [9], we only need to analyze a unit cell of a 2D periodic structure. The cross-section of a generic multilayer structure, which is periodic in the and directions, is displayed in Figure 1. We denote each homogeneous layer by , , where is the number of layers. , are the upper and lower half spaces, respectively. The interface between the and regions, denoted by , is located at . All 2D periodic thin conducting patterns are located at these interfaces (although some interfaces can be simple dielectric interfaces without any conducting patterns). Each homogeneous layer has a thickness of . The theory presented in this paper is applicable to two types of periodic structures:Type A: has a conducting pattern and/or a dielectric interface at ;Type B: has a PEC ground plane at .
An FSS with multiple printed metal patterns is an example of a Type A structure. A multilayer AMC surface is a Type B structure.
The following notations are adopted throughout this paper.(a): the incident electric field in Region .(b): the reflected electric field in Region if the top conducting pattern at were a complete PEC without any apertures. Note that, at , and .(c): the component of the equivalent magnetic currents at the top interface of region .(d): the component of the equivalent magnetic currents at the bottom interface of region .(e): the electric field created by and in Region .(f): the electric field in Region due to all equivalent magnetic currents.(g): the total electric field in Region .(h): is the component of the total magnetic vector potential in Region . and are the components of the magnetic vector potential created by the top and bottom equivalent magnetic currents, respectively, in Region .(i), : the permittivity and permeability of the dielectric in region .
The total field in region is given by where the scattered field is totally due to the magnetic currents at the bottom interface, since there is no top interface for region . In other regions, In Type B structures, .
In general, the conducting patterns on each interface have different periodicity from each other. In such cases, a larger unit cell with an effective periodicity is defined. The unit cell of a multilayer structure is considered as a group of coupled regions. Each region is bounded by Bloch-Floquet periodic boundaries on four sides. Let us take Region , , for example, as shown in Figure 2. The range of is . In Figure 2, without any loss of generality, we have assumed that and to simplify the explanation. The apertures in the top and bottom interfaces ( and ) of Region are replaced by equivalent magnetic currents on PECs. The top equivalent magnetic currents now appear in the aperture areas of Surface , just below the PEC at , while the bottom equivalent magnetic currents (, ) exist in the apertures of Surface , just above the bottom PEC at .
2.2. Expansion of the Coupled Fields in a Unit Cell of a Single-Layer Periodic Structure
In order to develop CFEM step by step, this subsection illustrates how to expand fields in each homogeneous region of a unit cell in a single-layer structure. For this purpose, we first transform the unit cell to two waveguides. The unit cell, shown in Figure 4(a), is separated into two regions, Region and Region , as shown in Figure 4(b). The metal pattern at is replaced with a PEC and equivalent magnetic currents in the aperture areas of both regions. Each region has Bloch-Floquet periodic boundaries on four sides.
As shown in Figure 4(b), the range of Region is . Its top is open and its bottom has equivalent magnetic currents, and , just above a PEC. In Type A structures, the range of Region is with an open boundary at the bottom. In Type B structures, the range of Region is with a PEC at the bottom as shown in Figure 4(b). In both cases, the top boundary of Region 1 is a PEC at . There exist equivalent magnetic currents, , , just below this PEC.
Next we apply image theory to Region and remove the PEC surface at , by doubling the magnetic currents. This transforms Region to the “waveguide” shown in Figure 3(a), which has Bloch-Floquet periodic boundaries on four sides, opened at two ends and excited by at . When image theory is applied to Region of a Type A structure, it is transformed to a similar waveguide, but excited by at . Region of a Type B structure is transformed to the waveguide shown in Figure 3(c), which is excited by at and their images at , . In our theory, we exploit the fact that the excitation in Figure 3(c) can be realized by superposing the excitation in Figure 3(c) at an infinite series of coordinates.
(a)
(b)
(c)
(d)
To find the electromagnetic fields produced by magnetic currents in Region and Region [10], we solve the following wave equations for the magnetic vector potential () in the two corresponding waveguides: where and represent the coordinate of each layer of equivalent magnetic currents (i.e., excitations in Figures 3(a) and 3(c)). is infinite at , otherwise it is zero.
The general solutions of (3) and (4), the inhomogeneous, second-order partial differential equations, are: Here, we select as eigenmodes of the waveguide bounded by Bloch-Floquet periodic boundary conditions. The sum of these modal functions, weighted with coefficients , gives the full-wave solution including all propagating and evanescent waves. The subscript represents the sequence of modal functions given below. These modal functions are the same in both regions: where , are the and components of the incident wave vector ; is the lattice constant along the direction; is the lattice constant along the direction; and are integers from to . In the special case of normal incidence, .
By substituting in the wave equation, one can prove that the in (7) is given by where For a Type A structure, is given by For a Type B structure, is given by To shorten the mathematical expressions, let us define two new functions, and , to represent the 3D field modes in (7) and their derivatives: where At , The electric and magnetic fields produced by the magnetic currents in each region are obtained by substituting (7) and (8) into (5) and (6):
2.3. Coupled-Field Expansion for a Unit Cell of a Multilayer Periodic Structure
Region , shown in Figure 2, is the generalized case of Region and Region in Section 2.2. This subsection extends the field expansion in Section 2.2 to the multilayer case. Let us first apply image theory to Region . Its top equivalent magnetic currents , are doubled when the top PEC at is removed. The same applies to its bottom equivalent magnetic currents , when the bottom PEC at is removed. In addition, multiple images of these currents appear in an infinite number of surfaces. This transforms Region to a waveguide with Bloch-Floquet periodic boundaries at the four sides, excited by and at and a series of their images at , as well as , at and a series of their images at , where . This excitation scheme for a typical Region is shown in Figure 3(d).
As the top and bottom regions are different from a typical region, shown in Figure 2, we need to consider four types of “waveguides”, as follows.(1)For , Region leads to the “waveguide” in Figure 3(a) because the top of this region is open and hence there are no magnetic currents on the top. This has only one layer of magnetic currents.(2)In Type A structures, for , Region leads to the “waveguide” in Figure 3(b) because the bottom of this region is open. In Type B structures, there is no field in this region as Surface is a PEC.(3)In Type B structures, for , Region leads to the “waveguide” in Figure 3(c) because Surface is a complete PEC without apertures and hence there are no equivalent magnetic currents on it.(4)All other typical regions have metallic patterns at both top and bottom, which leads to the generic “waveguide” in Figure 3(d).
Note that the last type of “waveguide” is required only for multilayer structures. We exploit the fact that the excitations in Figure 3(d) can be realized by shifting and superposing the excitations in Figures 3(a) and 3(b) at an infinite periodic series of coordinates.
To find the electromagnetic fields in Region produced by magnetic currents on Surfaces and , we solve the following wave equations for the magnetic vector potential () in the corresponding waveguides: where and , . is infinite at , otherwise, it is zero. is infinite at , otherwise it is zero. The series of and represents the coordinate of each layer of equivalent magnetic currents. (i.e., excitations in Figures 3(a), 3(b), 3(c), and 3(d)).
The magnetic potentials produced by the equivalent magnetic currents in the typical waveguide in Figure 3(d) can be expressed as a linear superposition of the potentials created by the equivalent magnetic currents at the two interfaces ( and ), and their corresponding series of images. So, let us separate them as Now, we expand these potentials as where and are defined in (9).
In (23)–(26), and are given by where To shorten the mathematical expressions, let us define four new functions, , , , and to represent the 3D field modes in (23) and (26), and their derivatives: where It can be shown that, at , Similarly, at ,
For further compactness, let us define three operators: where is an arbitrary function. , , and represent the integration over the aperture areas, metal patch areas, and the whole unit cell area at the interface , respectively. The operator satisfies the following orthogonality: Because the integration area of a unit cell is equal to the sum of its aperture area and patch area, In the following subsections, various interface continuities are enforced, in order to find the unknown coefficients and in (23)–(26).
Since Region is the generalized case of Region and in Section 2.2, we will enforce the boundary conditions on Surface .
2.4. Enforcing the Tangential Electric Field to Be Continuous across Apertures in Surface
Let us consider Region , and apply operator (35) to both sides of (48) at . The result is Due to the orthogonality in (36), only the , term is nonzero on the right-hand side. Then, considering that and as given in (32), we have Therefore, Similarly, we can derive the following expressions by applying the operator (35) to (48) for Region and applying (35) to (53) for Region and Region : Let us consider the interface in a typical region (). The component of the total electric field on this plane is given by . Then, as directed magnetic currents do not produce an directed electric field on the same plane (see (26)), , and, therefore, at this plane at . Similarly, we have for Region at the same plane at . Since the tangential components of the total electric field are zero on the metal areas of the interface, and they are continuous in the aperture areas of the interface, we find that In region , the component of the total electric field is the sum of the components of the incident field (), the reflected field (), and the scattered field. However, at the interface, and cancel out due to the PEC, and hence the above relationships are valid also for region at interface .
Let us substitute (44) into (40)-(41) and make a simple rearrangement which leads to the following relationship between the coefficients for the two coupled regions and : where is a new unknown coefficient, which will replace both and eventually.
Similarly, considering in (45) and (42)-(43), we have Note that the above two equations establish a simple relationship between the coefficients in Region and those in Region at . We have exploited this relationship to introduce a new set of unknown coefficients and , which are common to both regions. Note that the coefficients and are associated with the th interface whereas the former coefficients are associated with both the regions and interfaces.
Let us substitute (46) and (47) in (23)–(26). The electric and magnetic fields produced by the magnetic currents in each region are obtained by substituting these new expressions for the magnetic vector potentials in (20) and (21):
2.5. Enforcing the Tangential Electric Field to Be Zero on the Metal Patch Areas in Surface
As the tangential electric field is zero on the metal patch areas of the interface, and . Hence, we can replace in (40) by as we only consider at . In other words, the integration area in (40) can be reduced to the aperture area: Similarly, one may prove that Let us replace in (56) using and expand in (56) using (48). We obtain Note that (57) leads to the same equation (60), when the above steps are followed.
Similarly, substitution of (47) and (53) in (58)-(59) leads to where Note that the same integral (62) appears in both (60)-(61).
2.6. Enforcing the Tangential Magnetic Field to Be Continuous across the Aperture Areas in Surface
The total in region is where and are the components of the magnetic field of incident and reflected waves in region . At the interface , due to the PEC, and hence The total in a typical region, , is Now let us apply operator to (64): Next, we consider the interface and Region below it. So, we set in (65) and apply the same operator : Since the magnetic field is continuous across the aperture at interface , . So, we subtract (67) from (66). After simply rearranging the sequence of terms, we obtain Note that due to (37). Similarly the continuity of the tangential component of the magnetic field across the apertures at the other typical interfaces (, ) leads to
Let us repeat this process with . At interface , the relationship is For any other typical interface , ,
Next, we express all terms inside the operator in the right-hand side of (68)–(71) in terms of unknowns . Let us start from in (69). Replacing in (55) with and then applying the operator, we obtain
Due to the orthogonality in (36), only the , term is nonzero on the right-hand side. Hence, where
Next and are expanded using (61). Similarly, we apply operator to in Region . The result is where when , , otherwise .
We similarly expand the other terms in (68)–(71) and express the right-hand side of all these equations in terms of .
2.7. Building a Linear System of Equations
The linear system will be based on the expanded versions of (68)–(71).
2.7.1. Intermediate Interfaces
The matrix relationship developed in this subsection applies to all intermediate, typical surfaces . For Type A structures, they range from and, for Type B structures, .
Let us define three column vectors to represent unknowns in (69) and (71): where and , are from .
By substituting expansions (73)-(75), and so forth in (69) and (71), we obtain the linear subsystem for each intermediate layer: where The above expression describes the relationship between the unknowns for three adjacent surfaces. In this expression, , , , , , , , , , , , and are twelve submatrices. Their elements are where , and , , are given by
2.7.2. Last Interface ()
When the derivation in the above subsection is repeated for the last interface of a Type A structure, that is, , the result is This involves only two unknown vectors because there are no magnetic currents below this surface.
By manipulating the above equation, can be expressed in terms of : Let us now define then (83) becomes Note that now we have a relationship connecting unknowns for two adjacent interfaces.
In a Type structure, the last surface with magnetic currents is . The above derivation is valid for this case too. In this case, ; the relationship between and is similar to (85).
2.7.3. Second-Last Interface
As an example, let us now consider the next interface, . The objective here is to derive a relationship between and , by exploiting (78) and (85). Let us first substitute in (78) and rewrite it as Now, we substitute (85) to replace . The result is After rearrangement, we get the following relationship between and : where In a Type B structure, the second-last surface with magnetic currents is . The above derivation is valid for this case too.
2.7.4. Intermediate Interface
Now, we revisit a typical intermediate surface . As shown above for the example, it is possible to express the unknown coefficients () in terms of by repeating the same steps. The result is where
2.7.5. Top Interface
This interface is different due to the existence of the incident field in Region . After expanding (68) and (70), we obtain where is defined as Using (90), in (92) can be replaced by . The result is Note that the only unknown here is , which can be obtained by solving this matrix equation. Once is found, all other unknowns can be found, if necessary, step by step, following (90).
When calculating the reflection characteristics of multilayer structures, only is required. In fact, if all modes other than the , are evanescent, one needs only to determine the surface reflection coefficients.
2.8. Truncating the Infinite Series
It may appear that it is a large task to solve the linear system (92) because its unknowns are the coefficients of 2D modes for two orthogonal polarizations. For example, if we truncate the range of and as , then we have a linear system with unknowns. However, not all these unknowns contribute to the convergence equally. In fact, converged results can be obtained with a much smaller number of unknowns, as described here. For some geometries, only tens of unknowns are sufficient. There are two factors to be considered when determining the minimum set of unknowns.
First, the width of the spatial spectrum of electric magnetic fields is mainly determined by the singularities of the electromagnetic fields around the edges, tips, and reentrant corners of conducting patterns. These singularities can be described as in general, where is the distance to the singular centres, such as aperture corners, edges and tips, and is the singularity exponent (e.g., for a singularity along the edge of an infinite thin metal patch). Its spatial spectrum is given by , where is a constant for a given [11]. Take the unit cell, shown in Figure 6, for example. Since its singularity is mainly along the edges of the metal patch, the singularity exponent is equal to 0.5 (according to [12]) and according to [11]. This spectrum is effectively truncated when we truncate the modes by limiting the range of and . For this type of singularity in a mm unit cell, when the number of modes selected (i.e., ) is 7, its truncation error is less than . There should be at least 15 modes to have an error less than . Hence, the range of and is limited to for the calculations presented here. This would normally lead to unknowns, or altogether in both and directions.
However, not all the unknowns contribute to convergence equally, and hence most of them can be safely dropped without compromising accuracy. According to [13], the convergence of the linear system (94) mainly depends on the right-side vector . It can be shown that the (th) elements of and are proportional to , which is determined by the geometry of the conducting pattern. The of the periodic square patch, depicted in Figure 6, is shown in Figure 5. We note that along diagonal directions decay much faster than those along the axes, that is, or . It is significant only along the axes and around the origin. This makes the on-axis unknowns and dominant over off-axis unknowns. Hence, most of the off-axis (i.e., and ) can be ignored in the linear system without significantly affecting accuracy.
(a)
(b)
The modes that can be ignored are identified as follows. The right-hand-side vector is generated for all modes first. Based on the value of the largest element, a threshold is set and modes corresponding to elements below this threshold are ignored. The matrix in (94) is generated only after this truncation. Results obtained by applying this truncation scheme are given in the next section.
3. Validation and Results
The reflection and transmission coefficients of several planar periodic structures, for both TE and TM incidence, have been calculated to validate the accuracy of this method. These results are compared with HFSS results in this section. In the following derivation, we assume that all modes other than and are evanescent, which is the case with all structures we have analyzed so far. When a TE incident wave impinges on a periodic surface with incident angles (, ), the reflection and transmission coefficients ( and ) for the , mode are obtained as where are the three coordinate components of the unit vector , which is perpendicular to both the propagation direction vector and the surface normal direction : Similarly, for TM-polarized incident waves, the reflection and transmission coefficients are where For convenience and without any loss of generality, we chose to calculate the reflection and transmission coefficients at , . For the special case of normal incidence, For validation, we compare CFEM results with HFSS results because HFSS is readily available and is considered an industry standard software. In our HFSS model, the unit cell was built with Master/Slave boundaries on the four sides and a radiation boundary on the top (and also the bottom in the case of FSS) to truncate the computational volume. The model was excited by a TEM plane wave, instead of a wave port to define incident waves, which may impinge at oblique angles. In this case, the built-in calculation in HFSS cannot be used to calculate and . We instead used two other methods to calculate them. (i) Far-field method: we integrate the scattered field over a surface normal to the incident direction in the far field and find the average coefficient over the surface area. (ii) Near-field method: we first determine and by integrating the total field over the aperture area, according to and then calculate and using (95)–(97). In the first method, an average reflection coefficient is obtained directly using fields in the far-field region whereas, in the second method, the same is predicated using the aperture fields in the near-field region. We found that both far-field and near-field methods gave almost the same results.
In the following, we compare and calculated using CFEM, HFSS (far-field and near-field) methods. The first example is a Type B square patch AMC surface, whose geometry is sketched in Figure 6. The reflection phase of this square patch AMC, illuminated by a uniform plane wave incident in the normal direction (, ), has been calculated. As shown in Figure 7, excellent convergence was achieved with 593 Bloch-Floquet modes in the direction and the same number of Bloch-Floquet modes in the direction. The overall CPU time for the calculation of each frequency point was 40 seconds in a 3.0 GHz Pentium 4 processor. The CFEM results show good agreement with the HFSS results obtained using both near-field and far-field methods.
To assess the effect of the substrate thickness on the convergence, we reduced the substrate thickness in Figure 6 to 0.5 mm. Figure 8 shows that the method still converges with modes. This is because we analytically describe the effect of the substrate thickness in our model.
The TE and TM reflection phase for oblique incidence (, ) has been obtained with modes using CFEM. As shown in Figures 9 and 10, CFEM results again show good agreement with HFSS results. With the change of the incident angle from normal to oblique, the in-phase TE reflection frequency moves up slightly. For the angle of incidence , the AMC band for TM incidence is shifted upwards compared to TE incidence.
From Figures 7–10, we see there are slight discrepancies between the HFSS near-field and far-field results. These discrepancies are due to the slightly different numerical accuracy of the far-field and near-field values. These differences are not practically significant and hence both near-field and far-field methods are trustworthy.
The next example refers to a reflection-type circular polarizer, sketched in Figure 11, presented previously in [14]. When a wave linearly polarized in the diagonal direction is normally incident on this surface, a circularly polarized wave is reflected. The parameters of the unit cell are mm, mm, mm, mm, mm, mm, mm, mm; the thickness of the substrate is mm; the relative permittivity of the substrate is 4.4. An incident wave linearly polarized at to the axis has equal and components (, ). When all the modes other than , are evanescent, the reflection coefficients for and polarizations are This circular polarizer was designed such that the reflection coefficient magnitudes ( and ) are close to dB over a wide band (9.5–14.5 GHz) and the difference between the phases of the reflection coefficients () is close to 90° over the same band, as shown in Figure 12. The overall CPU time for the calculation of each frequency point was 180 seconds in a 3.0 GHz Pentium 4 processor and the memory required was 1.4 GBytes. The results obtained using CFEM are again in good agreement with HFSS results.
(a)
(b)
(c)
4. Conclusions
In this paper, a new, field-based, seminumerical method for analyzing single-layer and multilayer planar periodic structures has been proposed. Unlike many fully numerical methods, CFEM provides a physical insight as a result of the nature of its output, that is, field distribution in terms of modes. Software based on fully numerical methods (e.g., HFSS) provides values of the field at discrete sample points only.
Numerical results obtained using the CFEM for two 2D periodic structures agree well with HFSS results. One may note that, for AMC surfaces with thick (mm) and thin (mm) substrates, the method requires a similar number of unknowns to obtain convergent results. As demonstrated in Figures 7–10, the comparison between CFEM and HFSS results for both normal and oblique incidence is very good. Even for the complex geometry of the two-layer planar periodic metallic pattern (shown in Figure 12) that includes most types of singularities, accurate results have been obtained.
Acknowledgment
The authors would like to take this opportunity to thank ac3 and PETSC for their support with high-performance and scientific computations.