# A Coupled-Field Expansion Method for Single-Layer and Multilayer Planar Periodic Structures

Authors: Yading Li; Karu P. Esselle; Lvqian Zhang
Journal: International Journal of Antennas and Propagation
 (2012)

---

## 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), thex- 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.


---

## Body

## 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. Thex- and y-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 Section2.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 x- and y-directions, is displayed in Figure 1. We denote each homogeneous layer by Ri, i=1,…,n, where n is the number of layers. R0, Rn+1 are the upper and lower half spaces, respectively. The interface between the Ri-1 and Ri regions, denoted by Si, is located at z=zi. 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 Ri has a thickness of hi=zi+1-zi. 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 atz=zn+1;Type B: has a PEC ground plane atz=zn+1.Figure 1
Cross-section of a generic multilayer 2D periodic structure.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)
E⃗i=Exix̂+Eyiŷ+Eziẑ: the incident electric field in Region R0.(b)
E⃗r=Exrx̂+Eyrŷ+Ezrẑ: the reflected electric field in Region R0 if the top conducting pattern at z=z1 were a complete PEC without any apertures. Note that, at z=z1, Exi+Exr=0 and Eyi+Eyr=0.(c)
Mx(Ri,Si): the x-component of the equivalent magnetic currents at the top interface Si of region Ri.(d)
Mx(Ri,Si+1): the x-component of the equivalent magnetic currents at the bottom interface Si+1 of region Ri.(e)
E⃗(Ri)x=Ex(Ri)xx̂+Ey(Ri)xŷ+Ez(Ri)xẑ: the electric field created by Mx(Ri,Si) and Mx(Ri,Si+1) in Region Ri.(f)
E⃗(Ri)s=E⃗(Ri)x+E⃗(Ri)y: the electric field in Region Ri due to all equivalent magnetic currents.(g)
E⃗(Ri)tot=Ex(Ri)totx̂+Ey(Ri)totŷ+Ez(Ri)totẑ: the total electric field in Region Ri.(h)
𝒜(Ri)x=𝒜(Ri,Si)x+𝒜(Ri,Si+1)x: 𝒜(Ri)x is the x-component of the total magnetic vector potential in Region Ri. 𝒜(Ri,Si)x and 𝒜(Ri,Si+1)x are the x-components of the magnetic vector potential created by the top and bottom equivalent magnetic currents, respectively, in Region Ri.(i)
ϵ(Ri), μ(Ri): the permittivity and permeability of the dielectric in region Ri.The total field in regionR0 is given by(1)E⃗(Ri)tot=E⃗(Ri)s+E⃗i+E⃗r,
where the scattered field E⃗(Ri)s is totally due to the magnetic currents at the bottom interface, since there is no top interface for region R0. In other regions,(2)E⃗(Ri)tot=E⃗(Ri)s.
In Type B structures, E⃗(Ri+1)tot=0.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 RegionRi, (1≦i≦n), for example, as shown in Figure 2. The range of Ri is zi+1≦z≦zi. In Figure 2, without any loss of generality, we have assumed that zi+1=0 and zi=h to simplify the explanation. The apertures in the top and bottom interfaces (Si and Si+1) of Region Ri are replaced by equivalent magnetic currents on PECs. The top equivalent magnetic currents (Mx(Ri,Si),My(Ri,Si))  now appear in the aperture areas of Surface Si, just below the PEC at z=zi, while the bottom equivalent magnetic currents (Mx(Ri,Si+1), My(Ri,Si+1)) exist in the apertures of Surface Si+1, just above the bottom PEC at z=zi+1.Figure 2
Typical regionRi. It is assumed that zi+1=0, zi=h.


### 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 Figure4(a), is separated into two regions, Region R0 and Region R1, as shown in Figure 4(b). The metal pattern at z=h 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 Figure4(b), the range of Region R0 is z∈[h,+∞). Its top is open and its bottom has equivalent magnetic currents, Mx(R0) and My(R0), just above a PEC. In Type A structures, the range of Region R1 is z∈(-∞,h] with an open boundary at the bottom. In Type B structures, the range of Region R1 is z∈[0,h] 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 z=h. There exist equivalent magnetic currents, Mx(R1), My(R1), just below this PEC.Next we apply image theory to RegionR0 and remove the PEC surface at z=h, by doubling the magnetic currents. This transforms Region R0 to the “waveguide” shown in Figure 3(a), which has Bloch-Floquet periodic boundaries on four sides, opened at two ends and excited by 2M⃗(R0)=2Mx(R0)x̂+2M⃗y(R0)ŷ at z=h. When image theory is applied to Region R1 of a Type A structure, it is transformed to a similar waveguide, but excited by 2M⃗(R1)=2Mx(R1)x̂+2M⃗y(R1)ŷ at z=h. Region R1 of a Type B structure is transformed to the waveguide shown in Figure 3(c), which is excited by 2M⃗(R1) at z=h and their images at z=(2p+1)h, p=…,-2,-1,1,2,…. 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 z coordinates.Equivalent “waveguides” with Bloch periodic boundaries at four sides, opened at two ends; (a) excited by equivalent magnetic currents atz=h; (b) excited by equivalent magnetic currents at z=0; (c) excited by equivalent magnetic currents and their images at  z=(2p+1)h, p=…,-1,0,1,…; (d) excited by equivalent magnetic currents and their images at z=(2p+1)h and z=(2p)h, where p=…,-1,0,1,….
(a)(b)(c)(d)Figure 4
(a) A unit cell. (b) Separation of a Type B unit cell into two regions. RegionR0 is open at the top and has a PEC at the bottom. Equivalent magnetic currents (shown in light color at z=h) exist just above the bottom PEC. Region R1 has a PEC surface at the top and equivalent magnetic currents (shown in light color at z=h) exist just below the top PEC. For Type A structures, the bottom of Region R1 is open (not shown here). For Type B structures, it is a PEC (as shown).To find the electromagnetic fields produced by magnetic currents in RegionR0 and Region R1 [10], we solve the following wave equations for the magnetic vector potential (𝒜) in the two corresponding waveguides:(3)∇2Ax+k2Ax=-2Mxjωμδ(z-zp+1),(4)∇2Ay+k2Ay=-2Myjωμδ(z-zp+1),(5)E⃗s=-jωμ∇×A⃗,(6)H⃗s=∇×∇×A⃗,
where 𝒜⃗=𝒜xx̂+𝒜yŷ and zp+1=(2p+1)h represent the z coordinate of each layer of equivalent magnetic currents (i.e., excitations in Figures 3(a) and 3(c)). δ(z-zp+1) is infinite at z=zp+1, otherwise it is zero.The general solutions of (3) and (4), the inhomogeneous, second-order partial differential equations, are:(7)A(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)xFX(m)FY(n)FZ(R0,R1)(m,n),(8)A(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)yFX(m)FY(n)FZ(R0,R1)(m,n).
Here, we select FX(m)FY(n) as eigenmodes of the waveguide bounded by Bloch-Floquet periodic boundary conditions. The sum of these modal functions, weighted with coefficients amn(R0,R1)x,y, gives the full-wave solution including all propagating and evanescent waves. The subscript mn represents the sequence of modal functions given below. These modal functions are the same in both regions:(9)FX(m)=e-jkmx,km=2mπa+kx,FY(n)=e-jkny,kn=2nπb+ky,
where kx, ky are the x and y components of the incident wave vector k⃗(R0); a is the lattice constant along the x direction;  b is the lattice constant along the y direction; m and n are integers from -∞ to +∞. In the special case of normal incidence, kx=ky=0.By substituting in the wave equation, one can prove that theFZ(R0)(m,n) in (7) is given by(10)FZ(R0)(m,n)=e-γmn(R0)(z-h),
where(11)(γmn(R0,R1))2=km2+kn2-(k(R0,R1))2,(k(R0,R1))2=ω2ϵ0ϵr(R0,R1)μ0μr(R0,R1).
For a Type A structure, FZ(R1)(m,n) is given by(12)FZ(R1)(m,n)=eγmn(R1)(z-h).
For a Type B structure, FZ(R1)(m,n) is given by(13)FZ(R1)(m,n)=∑p=0+∞[eγmn(R1)(z-(2p+1)h)+e-γmn(R1)(z+(2p+1)h)].
To shorten the mathematical expressions, let us define two new functions, F(R0,R1) and Fd(R0,R1), to represent the 3D field modes in (7) and their derivatives:(14)F(R0,R1)(m,n)=FX(m)FY(n)FZ(R0,R1)(m,n),Fd(R0,R1)(m,n)=FX(m)FY(n)FZd(R0,R1)(m,n),
where(15)FZd(R0,R1)(m,n)=-∂FZ(R0,R1)(m,n)γmn(R0,R1)∂z.
At z=h,(16)FZd(R0)=1,FZd(R1)=-1.
The electric and magnetic fields produced by the magnetic currents in each region are obtained by substituting (7) and (8) into (5) and (6):(17)Ex(R0,R1)x=0,Ey(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)x(jωμ(R0,R1)γmn(R0,R1))Fd(R0,R1)(m,n),Hx(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)x((k(R0,R1))2-(km)2)×F(R0,R1)(m,n),Hy(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)x(-kmkn)F(R0,R1)(m,n),Ex(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)y(-jωμ(R0,R1)γmn(R0,R1))×Fd(R0,R1)(m,n),Ey(R0,R1)y=0,Hx(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)y(-kmkn)F(R0,R1)(m,n),Hy(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)y((k(R0,R1))2-(kn)2)×F(R0,R1)(m,n).


### 2.3. Coupled-Field Expansion for a Unit Cell of a Multilayer Periodic Structure

RegionRi, shown in Figure 2, is the generalized case of Region R0 and Region R1 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 Ri. Its top equivalent magnetic currents Mx(Ri,Si), My(Ri,Si) are doubled when the top PEC at z=zi is removed. The same applies to its bottom equivalent magnetic currents Mx(Ri,Si+1), My(Ri,Si+1) when the bottom PEC at z=zi+1 is removed. In addition, multiple images of these currents appear in an infinite number of surfaces. This transforms Region Ri to a waveguide with Bloch-Floquet periodic boundaries at the four sides, excited by 2Mx(Ri,Si) and 2My(Ri,Si) at z=zi and a series of their images at z=zi+(2p)hi, as well as 2Mx(Ri,Si+1), 2My(Ri,Si+1) at z=zi+1 and a series of their images at z=zi+1+(2p)hi, where p=…,-2,-1,1,2,…. This excitation scheme for a typical Region i is shown in Figure 3(d).As the top and bottom regions are different from a typical region, shown in Figure2, we need to consider four types of “waveguides”, as follows.(1)
Fori=0, Region R0 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, fori=n+1, Region Rn+1 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 Sn+1 is a PEC.(3)
In Type B structures, fori=n, Region Rn leads to the “waveguide” in Figure 3(c) because Surface Sn+1 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 Figure3(d).Note that the last type of “waveguide” is required only for multilayer structures. We exploit the fact that the excitations in Figure3(d) can be realized by shifting and superposing the excitations in Figures 3(a) and 3(b) at an infinite periodic series of z-coordinates.To find the electromagnetic fields in RegionRi produced by magnetic currents on Surfaces Si and Si+1, we solve the following wave equations for the magnetic vector potential (𝒜(Ri)) in the corresponding waveguides:(18)∇2A(Ri)x+(k(Ri))2A(Ri)x=-2Mx(Ri,Si)δ(z-zi-zp)+2Mx(Ri,Si+1)δ(z-zi+1-zp)jωμ(Ri),(19)∇2A(Ri)y+(k(Ri))2A(Ri)y=-2My(Ri,Si)δ(z-zi-zp)+2My(Ri,Si+1)δ(z-zi+1-zp)jωμ(Ri),(20)E⃗(Ri)s=-jωμ(Ri)∇×A⃗(Ri),(21)H⃗(Ri)s=∇×∇×A⃗(Ri),
where 𝒜⃗(Ri)=𝒜(Ri)xx̂+𝒜(Ri)yŷ and zp=(2p)hi, p=…,-2,-1,0,1,2,…. δ(z-zi-(2p)hi) is infinite at z=zi+(2p)hi, otherwise, it is zero. δ(z-zi+1-(2p)hi) is infinite at z=zi+1+(2p)hi, otherwise it is zero. The series of z=zi+zp and z=zi+1+zp   represents the z 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 Figure3(d) can be expressed as a linear superposition of the potentials created by the equivalent magnetic currents at the two interfaces (z=zi and z=zi+1), and their corresponding series of images. So, let us separate them as(22)A(Ri)y=A(Ri,Si+1)y+A(Ri,Si)y,A(Ri)x=A(Ri,Si+1)x+A(Ri,Si)x.
Now, we expand these potentials as(23)A(Ri,Si+1)y=∑m,n=-∞+∞amn(Ri,Si+1)yFX(m)FY(n)FZ(Ri,Si+1)(m,n),(24)A(Ri,Si)y=∑m,n=-∞+∞amn(Ri,Si)yFX(m)FY(n)FZ(Ri,Si)(m,n),(25)A(Ri,Si+1)x=∑m,n=-∞+∞amn(Ri,Si+1)xFX(m)FY(n)FZ(Ri,Si+1)(m,n),(26)A(Ri,Si)x=∑m,n=-∞+∞amn(Ri,Si)xFX(m)FY(n)FZ(Ri,Si)(m,n),
where FX(m) and FY(n) are defined in (9).In (23)–(26), FZ(Ri,Si+1)(m,n) and  FZ(Ri,Si)(m,n) are given by(27)FZ(Ri,Si+1)(m,n)=∑p=0+∞[eγmn(Ri)(z-zi+1-(2p)hi)+e-γmn(Ri)(z-zi+1+(2p)hi)],FZ(Ri,Si)(m,n)=∑p=0+∞[eγmn(Ri)(z-zi-(2p)hi)+e-γmn(Ri)(z-zi+(2p)hi)],
where(28)(γm,n(Ri))2=(km)2+(kn)2-(k(Ri))2,(k(Ri))2=ω2μ(Ri)ϵ(Ri).
To shorten the mathematical expressions, let us define four new functions, F(Ri,Si+1), F(Ri,Si), Fd(Ri,Si+1),  and Fd(Ri,Si)  to represent the 3D field modes in (23) and (26), and their derivatives:(29)F(Ri,Si+1)(m,n)=FX(m)FY(n)FZ(Ri,Si+1)(m,n),F(Ri,Si)(m,n)=FX(m)FY(n)FZ(Ri,Si)(m,n),Fd(Ri,Si+1)(m,n)=FX(m)FY(n)FZd(Ri,Si+1)(m,n),Fd(Ri,Si)(m,n)=FX(m)FY(n)FZd(Ri,Si)(m,n),
where(30)FZd(Ri,Si+1)(m,n)=-1γmn(Ri)⋅∂FZ(Ri,Si+1)(m,n)∂z,FZd(Ri,Si)(m,n)=1γmn(Ri)⋅∂FZ(Ri,Si)(m,n)∂z.
It can be shown that, at z=zi,(31)FZd(Ri,Si+1)=0,FZd(Ri,Si)=1.
Similarly, at z=zi+1,(32)FZd(Ri,Si+1)=1,FZd(Ri,Si)=0.For further compactness, let us define three operators:(33)Γa(Si)m,n[f(x)]=∬apf(x)FX*(m)FY*(n)dxdy,(34)Γp(Si)m,n[f(x)]=∬patchf(x)FX*(m)FY*(n)dxdy,(35)Γu(Si)m,n[f(x)]=∬unitf(x)FX*(m)FY*(n)dxdy,
where f(x) is an arbitrary function. ∬ap, ∬patch, and ∬unit  represent the integration over the aperture areas, metal patch areas, and the whole unit cell area at the interface Si, respectively. The operator Γu(Si)m,n satisfies the following orthogonality:(36)Γu(Si)m,n[FX(m′)FY(n′)]={ab,m′=m,n′=n,0,otherwise.
Because the integration area of a unit cell is equal to the sum of its aperture area and patch area,(37)Γu(Si)m,n[f(x)]=Γa(Si)m,n[f(x)]+Γp(Si)m,n[f(x)].
In the following subsections, various interface continuities are enforced, in order to find the unknown coefficients amn(Ri,Si+1)x,y and amn(Ri,Si)x,y in (23)–(26).Since RegionRi is the generalized case of Region R0 and R1 in Section 2.2, we will enforce the boundary conditions on Surface Si.


### 2.4. Enforcing the Tangential Electric Field to Be Continuous across Apertures in SurfaceSi

Let us consider RegionRi-1, and apply operator (35) to both sides of (48) at z=zi. The result is(38)Γu(Si)m,n[Ex(Ri-1)y]=Γu(Si)m,n[∑m′,n′=-∞+∞(jωμ(Ri-1)γm′n′(Ri-1))⋅(am′n′(Ri-1,Si-1)yFd(Ri-1,Si-1)-am′n′(Ri-1,Si)yFd(Ri-1,Si))∑m′,n′=-∞+∞].
Due to the orthogonality in (36), only the m′=m, n′=n term is nonzero on the right-hand side. Then, considering that Fd(Ri-1,Si-1)|z=zi≡0 and Fd(Ri-1,Si)|z=zi≡1 as given in (32), we have(39)Γu(Si)m,n[Ex(Ri-1)y]=amn(Ri-1,Si)y(-jωμ(Ri-1)γmn(Ri-1))ab.
Therefore,(40)amn(Ri-1,Si)y=Γu(Si)m,n[Ex(Ri-1)y]-jωμ(Ri-1)γmn(Ri-1)ab.
Similarly, we can derive the following expressions by applying the operator (35) to (48) for Region Ri and applying (35) to (53) for Region Ri-1 and Region Ri:
(41)amn(Ri,Si)y=Γu(Si)m,n[Ex(Ri)y]jωμ(Ri)γmn(Ri)ab,(42)amn(Ri-1,Si)x=Γu(Si)m,n[Ey(Ri-1)x]jωμ(Ri-1)γmn(Ri-1)ab,(43)amn(Ri,Si)x=Γu(Si)m,n[Ey(Ri)x]-jωμ(Ri)γmn(Ri)ab.
Let us consider the interface Si in a typical region Ri-1 (i-1≠0,n). The x-component of the total electric field on this plane is given by Ex(Ri-1)tot=Ex(Ri-1)x+Ex(Ri-1)y. Then, as x-directed magnetic currents do not produce an x-directed electric field on the same plane (see (26)), Ex(Ri-1)x=0, and, therefore, Ex(Ri-1)tot=Ex(Ri-1)y at this plane at z=zi. Similarly, we have Ex(Ri)tot=Ex(Ri)y for Region Ri at the same plane at z=zi. 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(44)Γu(Si)m,n[Ex(Ri-1)y]=Γu(Si)m,n[Ex(Ri)y],(45)Γu(Si)m,n[Ey(Ri-1)x]=Γu(Si)m,n[Ey(Ri)x].
In region R0, the x-component of the total electric field is the sum of the x-components of the incident field (Exi), the reflected field (Exr), and the scattered field. However, at the S1 interface, Exi and Exr cancel out due to the PEC, and hence the above relationships are valid also for region R0 at interface S1.Let us substitute (44) into (40)-(41) and make a simple rearrangement which leads to the following relationship between the amn coefficients for the two coupled regions Ri and Ri+1:(46)bmn(Si)y=amn(Ri-1,Si)y(-jωμ(Ri-1)γmn(Ri-1))=amn(Ri,Si)y(jωμ(Ri)γmn(Ri)),
where bmn(Si)y is a new unknown coefficient, which will replace both amn(Ri) and amn(Ri-1) eventually.Similarly, consideringEy in (45) and (42)-(43), we have(47)bmn(Si)x=amn(Ri-1,Si)x(jωμ(Ri-1)γmn(Ri-1))=amn(Ri,Si)x(-jωμ(Ri)γmn(Ri)).
Note that the above two equations establish a simple relationship between the coefficients in Region Ri-1 and those in Region Ri at z=zi. We have exploited this relationship to introduce a new set of unknown coefficients bmn(Si)y and bmn(Si)x, which are common to both regions. Note that the coefficients bmn(Si)x and bmn(Si)y are associated with the ith interface whereas the former coefficients amn(Ri,Sj) 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):(48)Ex(Ri)y=∑m′,n′=-∞+∞bm′n′(Si+1)yFd(Ri,Si+1)(m′,n′)+∑m′,n′=-∞+∞bm′n′(Si)yFd(Ri,Si)(m′,n′),(49)Ey(Ri)y=0,(50)Hx(Ri)y=∑m′,n′=-∞+∞bm′n′(Si+1)y(-km′kn′)F(Ri,Si+1)(m′,n′)-jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)y(-km′kn′)F(Ri,Si)(m′,n′)jωμ(Ri)γm′n′(Ri),(51)Hy(Ri)y=∑m′,n′=-∞+∞bm′n′(Si+1)y((k(Ri))2-kn′2)F(Ri,Si+1)(m′,n′)-jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)y((k(Ri))2-kn′2)F(Ri,Si)(m′,n′)jωμ(Ri)γm′n′(Ri),(52)Ex(Ri)x=0,(53)Ey(Ri)x=∑m′,n′=-∞+∞bm′n′(Si+1)xFd(Ri,Si+1)(m′,n′)+∑m′,n′=-∞+∞bm′n′(Si)xFd(Ri,Si)(m′,n′),(54)Hx(Ri)x=∑m′,n′=-∞+∞bm′n′(Si+1)x((k(Ri))2-km′2)F(Ri,Si+1)(m′,n′)jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)x((k(Ri))2-km′2)F(Ri,Si)(m′,n′)-jωμ(Ri)γm′n′(Ri),(55)Hy(Ri)x=∑m′,n′=-∞+∞bm′n′(Si+1)x(-km′kn′)F(Ri,Si+1)(m′,n′)jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)x(-km′kn′)F(Ri,Si)(m′,n′)-jωμ(Ri)γm′n′(Ri).


### 2.5. Enforcing the Tangential Electric Field to Be Zero on the Metal Patch Areas in SurfaceSi

As the tangential electric field is zero on the metal patch areas of the interface,Γp(Si)m,n[Ex(Ri-1)y]=0 and Γu(Si)m,n[Ex(Ri-1)y]=Γa(Si)m,n[Ex(Ri-1)y]. Hence, we can replace Γu(Si)m,n in (40) by Γa(Si)m,n as we only consider Ex(Ri-1)y at z=zi. In other words, the integration area in (40) can be reduced to the aperture area:(56)amn(Ri-1,Si)y=Γa(Si)m,n[Ex(Ri-1,Si)y]-jωμ(Ri-1)γmn(Ri-1)ab.
Similarly, one may prove that(57)amn(Ri,Si)y=Γa(Si)m,n[Ex(Ri)y]jωμ(Ri)γmn(Ri)ab,(58)amn(Ri-1,Si)x=Γa(Si)m,n[Ey(Ri-1)x]jωμ(Ri-1)γmn(Ri-1)ab,(59)amn(Ri,Si)x=Γa(Si)m,n[Ey(Ri)x]-jωμ(Ri)γmn(Ri)ab.
Let us  replace amn(Ri-1,Si)y(-jωμ(Ri-1)γmn(Ri-1)) in (56) using bmn(Si)y and expand Ex(Ri-1,Si)y in (56) using (48). We obtain(60)bmn(Si)y=∑m′,n′=-∞+∞bm′n′(Si)yΛap(Si)(m,n,m′,n′)ab.
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(61)bmn(Si)x=∑m′,n′=-∞+∞bm′n′(Si)xΛap(Si)(m,n,m′,n′)ab,
where(62)Λap(Si)(m,n,m′,n′)=Γa(Si)m,n[FX(m′)FY(n′)].
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 SurfaceSi

The totalHy in region R0 is(63)Hy(R0)tot=Hy(R0)x+Hy(R0)y+Hyi+Hyr,
where Hyi and Hyr are the y-components of the magnetic field of incident and reflected waves in region R0. At the interface S1, Hyr=Hyi due to the PEC, and hence(64)Hy(R0)tot=Hy(R0)x+Hy(R0)y+2Hyi.
The total Hy in a typical region, Ri, 1<i<n+1 is(65)Hy(Ri)tot=Hy(Ri)x+Hy(Ri)y.
Now let us apply operator Γa(S1)m,n to (64):(66)Γa(S1)m,n[Hy(R0)tot]=Γa(S1)m,n[Hy(R0)x+Hy(R0)y+2Hyi].
Next, we consider the interface S1 and Region R1 below it. So, we set i=1 in (65) and apply the same operator Γa(S1)m,n:(67)Γa(S1)m,n[Hy(R1)tot]=Γa(S1)m,n[Hy(R1)x+Hy(R1)y].
Since the magnetic field is continuous across the aperture at interface S1, Γa(S1)m,n[Hy(R0)tot]=Γa(S1)m,n[Hy(R1)tot]. So, we subtract (67) from (66). After simply rearranging the sequence of terms, we obtain(68)Γa(S1)m,n[-2Hyi]=Γa(S1)m,n[Hy(R0)x+Hy(R0)y-Hy(R1)x-Hy(R1)y]=(Γu(S1)m,n-Γp(S1)m,n)[Hy(R0)x+Hy(R0)y-Hy(R1)x-Hy(R1)y].
Note that Γa(S1)m,n=Γu(S1)m,n-Γp(S1)m,n due to (37). Similarly the continuity of the tangential component of the magnetic field across the apertures at the other typical interfaces (Si, i>1) leads to(69)0=Γa(Si)m,n[Hy(Ri-1)x+Hy(Ri-1)y-Hy(Ri)x-Hy(Ri)y]=(Γu(Si)m,n-Γp(Si)m,n)×[Hy(Ri-1)x+Hy(Ri-1)y-Hy(Ri)x-Hy(Ri)y].Let us  repeat this process withHx. At interface S1, the relationship is(70)Γa(S1)m,n[-2Hxi]=Γa(S1)m,n[Hx(R0)x+Hx(R0)y-Hx(R1)x-Hx(R1)y]=(Γu(S1)m,n-Γp(S1)m,n)[Hx(R0)x+Hx(R0)y-Hx(R1)x-Hx(R1)y].
For any other typical interface Si, i>1,(71)0=Γa(Si)m,n[Hx(Ri-1)x+Hx(Ri-1)y-Hx(Ri)x-Hx(Ri)y]=(Γu(Si)m,n-Γp(Si)m,n)×[Hx(Ri-1)x+Hx(Ri-1)y-Hx(Ri)x-Hx(Ri)y].Next, we express all terms inside theΓu(Si)m,n-Γp(Si)m,n operator in the right-hand side of (68)–(71) in terms of unknowns bmnx,y. Let us start from Γu(Si)m,n[Hy(Ri-1)x] in (69). Replacing i in (55) with i-1 and then applying the operator, we obtain(72)Γu(Si)m,n[Hy(Ri-1)x]=Γu(Si)m,n[∑m′,n′=-∞+∞bm′n′(Si)x(-km′kn′)F(Ri-1,Si)(m′,n′)jωμ(Ri-1)γm′n′(Ri-1)+∑m′,n′=-∞+∞bm′n′(Si-1)x(-km′kn′)F(Ri-1,Si-1)(m′,n′)-jωμ(Ri-1)γm′n′(Ri-1)].Due to the orthogonality in (36), only the m′=m, n′=n term is nonzero on the right-hand side. Hence,(73)Γu(Si)m,n[Hy(Ri-1)x]=(-kmkn)bmn(Si)xζ(Ri-1,Si)m,nabjωμ(Ri-1)γmn(Ri-1)+(-kmkn)bmn(Si-1)xζ(Ri-1,Si-1)m,nab-jωμ(Ri-1)γmn(Ri-1),
where(74)ζ(Ri-1,Si)m,n=FZ(Ri-1,Si)(m,n)|z=zi,ζ(Ri,Si)m,n=FZ(Ri,Si)(m,n)|z=zi.Nextbmn(Si)x and bmn(Si-1)x are expanded using (61). Similarly, we apply operator Γp(Si)m,n to Hy(Ri)x in Region Ri-1. The result is(75)Γp(Si)m,n[Hy(Ri-1)x]=∑m′,n′=-∞+∞(-km′kn′)bm′n′(Si)xζ(Ri-1,Si)m′,n′(Λ-Λap(Si)(m,n,m′,n′))jωμ(Ri-1)γm′n′(Ri-1)+∑m′,n′=-∞+∞(-km′kn′)bm′n′(Si-1)xζ(Ri-1,Si-1)m′,n′(Λ-Λap(Si)(m,n,m′,n′))-jωμ(Ri-1)γm′n′(Ri-1),
where Λ=ab when m=m′, n=n′, otherwise Λ=0.We similarly expand the other terms in (68)–(71) and express the right-hand side of all these equations in terms of bmnx,y.


### 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 InterfacesSi

The matrix relationship developed in this subsection applies to all intermediate, typical surfacesSi. For Type A structures, they range from i=2,…,n and, for Type B structures, i=2,…,n-1.Let us  define three column vectors to represent unknowns in (69) and (71):(76)B(Si-1)=[B(Si-1)yB(Si-1)x],B(Si)=[B(Si)yB(Si)x],B(Si+1)=[B(Si+1)yB(Si+1)x],
where(77)B(Si-1)y=[…,bm′n′(Si-1)y,…]′;B(Si-1)x=[…,bm′n′(Si-1)x,…]′,B(Si)y=[…,bm′n′(Si)y,…]′;B(Si)x=[…,bm′n′(Si)x,…]′,B(Si+1)y=[…,bm′n′(Si+1)y,…]′;B(Si+1)x=[…,bm′n′(Si+1)x,…]′,
and m′, n′ are from -∞→∞.By substituting expansions (73)-(75), and so forth in (69) and (71), we obtain the linear subsystem for each intermediate layer:(78)0=[ZSi(Si-1)ZSi(Si)ZSi(Si+1)]⋅[B(Si-1)B(Si)B(Si+1)],
where(79)ZSi(Si-1)=[OSi(Si-1)yPSi(Si-1)xQSi(Si-1)yLSi(Si-1)x],ZSi(Si)=[OSi(Si)yPSi(Si)xQSi(Si)yLSi(Si)x],ZSi(Si+1)=[OSi(Si+1)yPSi(Si+1)xQSi(Si+1)yLSi(Si+1)x].
The above expression describes the relationship between the unknowns for three adjacent surfaces. In this expression, OSi(Si-1)y, OSi(Si)y, OSi(Si+1)y, PSi(Si-1)x, PSi(Si)x, PSi(Si+1)x, QSi(Si-1)y, QSi(Si)y, QSi(Si+1)y, LSi(Si-1)x, LSi(Si)x, and LSi(Si+1)x are twelve submatrices. Their elements are(80)OSi(Si-1)y=τ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′-τ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,PSi(Si-1)x=-ξ(Ri-1)m,mζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′+ξ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,QSi(Si-1)y=ξ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,LSi(Si-1)x=-φ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′+φ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,OSi(Si)y=(τ(Ri)m,nζ(Ri,Si)m,n-τ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′-(τ(Ri)m′,n′ζ(Ri,Si)m′,n′-τ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,PSi(Si)x=-(ξ(Ri)m,nζ(Ri,Si)m,n-ξ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′+(ξ(Ri)m′,n′ζ(Ri,Si)m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,QSi(Si)y=(ξ(Ri)m,nζ(Ri,Si)m,n-ξ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′-(ξ(Ri)m′,n′ζ(Ri,Si)m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,LSi(Si)x=-(φ(Ri)m,nζ(Ri,Si)m,n-φ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′+(φ(Ri)m′,n′ζ(Ri,Si)m′,n′-φ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,OSi(Si+1)y=-τ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′+τ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,PSi(Si+1)x=ξ(Ri)m,mζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′-ξ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,QSi(Si+1)y=-ξ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′+ξ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,LSi(Si+1)x=φ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′-φ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,
where Λp(Si)m,n,m′,n′=Λ-Λap(Si)m,n,m′,n′, and ξ(Ri)m,n, τ(Ri)m,n, φ(Ri)m,n are given by(81)ξ(Ri)m,n=(-kmkn)jωμ(Ri)γmn(Ri),τ(Ri)m,n=((k(Ri))2-kn2)jωμ(Ri)γm,n(Ri),φ(Ri)m,n=((k(Ri))2-km2)jωμ(Ri)γm,n(Ri).


#### 2.7.2. Last Interface (Sn+1)

When the derivation in the above subsection is repeated for the last interfaceSn+1 of a Type A structure, that is, Sn+1, the result is(82)0=[ZSn+1(Sn)ZSn+1(Sn+1)]⋅[B(Sn)B(Sn+1)].
This involves only two unknown vectors because there are no magnetic currents below this surface.By manipulating the above equation,B(Sn+1) can be expressed in terms of B(Sn):(83)B(Sn+1)=[ZSn+1(Sn+1)]-1[-ZSn+1(Sn)B(Sn)].
Let us now define(84)Φn+1=-[ZSn+1(Sn+1)]-1ZSn+1(Sn),
then (83) becomes(85)B(Sn+1)=Φn+1B(Sn).
Note that now we have a relationship connecting unknowns for two adjacent interfaces.In a TypeB structure, the last surface with magnetic currents is Sn. The above derivation is valid for this case too. In this case, B(Sn+1)=0; the relationship between B(Sn) and B(Sn-1) is similar to (85).


#### 2.7.3. Second-Last InterfaceSn

As an example, let us now consider the next interface,Sn. The objective here is to derive a relationship between B(Sn) and B(Sn-1), by exploiting (78) and (85). Let us first substitute i=n in (78) and rewrite it as(86)ZSn(Sn)B(Sn)+ZSn(Sn+1)B(Sn+1)=-ZSn(Sn-1)B(Sn-1).
Now, we substitute (85) to replace B(Sn+1). The result is(87)[ZSn(Sn)+ZSn(Sn+1)Φn+1]B(Sn)=-ZSn(Sn-1)B(Sn-1).
After rearrangement, we get the following relationship between B(Sn) and B(Sn-1):(88)B(Sn)=ΦnB(Sn-1),
where(89)Φn=-[ZSn(Sn)+ZSn(Sn+1)Φn+1]-1⋅ZSn(Sn-1).
In a Type B structure, the second-last surface with magnetic currents is Sn-1. The above derivation is valid for this case too.


#### 2.7.4. Intermediate InterfaceSi

Now, we revisit a typical intermediate surfaceSi. As shown above for the i=n example, it is possible to express the unknown coefficients B(Si) (i=2,…,n+1) in terms of B(Si-1) by repeating the same steps. The result is(90)B(Si)=ΦiB(Si-1),
where(91)Φi=-[ZSi(Si)+ZSi(Si+1)Φi+1]-1⋅ZSi(Si-1).


#### 2.7.5. Top InterfaceS1

This interface is different due to the existence of the incident field in RegionR0. After expanding (68) and (70), we obtain(92)H(S1)=[ZS1(S1)ZS1(S2)]⋅[B(S1)B(S2)],
where H(S1) is defined as(93)H(S1)=[Hy(S1)iHx(S1)i],Hy(S1)i=[…,Γam,n[-2Hyi],…]′,Hx(S1)i=[…,Γam,n[-2Hxi],…]′.
Using (90), B(S2) in (92) can be replaced by B(S1). The result is(94)H(S1)=[ZS1(S1)+ZS1(S2)Φ2]⋅[B(S1)].
Note that the only unknown here is B(S1), which can be obtained by solving this matrix equation. Once B(S1) is found, all other unknowns can be found, if necessary, step by step, following (90).When calculating the reflection characteristics of multilayer structures, onlyB(S1) is required. In fact, if all modes other than the m=0, n=0 are evanescent, one needs only b00(S1) 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 m and n as [-15,…,0,…,15], then we have a linear system with 2×31×31≈2000 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 as1/R1-ν in general, where R is the distance to the singular centres, such as aperture corners, edges and tips, and ν is the singularity exponent (e.g., ν=0.5 for a singularity along the edge of an infinite thin metal patch). Its spatial spectrum is given by 2Ψ(ν)(k)-νcos(νπ/2), 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 Ψ(ν)=1.7725 according to [11]. This spectrum is effectively truncated when we truncate the modes by limiting the range of m and n. For this type of singularity in a 10mm unit cell, when the number of modes selected (i.e., maxm,n) is 7, its truncation error is less than 2.4%. There should be at least 15 modes to have an error less than 1%. Hence, the range of m and n is limited to [-15,…,0,…,15] for the calculations presented here. This would normally lead to 31×31 unknowns, or altogether 2×31×31 in both x- and y-directions.However, not all the2×31×31 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 [Hyi,Hxi]′. It can be shown that the ((m,n)th) elements of Hyi and Hxi are proportional to Λapm,n,0,0, which is determined by the geometry of the conducting pattern. The |Λapm,n,0,0| of the periodic square patch, depicted in Figure 6, is shown in Figure 5. We note that |Λapm,n,0,0| along diagonal directions decay much faster than those along the axes, that is, |Λapm,0,0,0| or |Λap0,n,0,0|. It is significant only along the axes and around the origin. This makes the on-axis unknowns bm0 and b0n dominant over off-axis unknowns. Hence, most of the off-axis bmn (i.e., m≠0 and n≠0) can be ignored in the linear system without significantly affecting accuracy.Figure 5
|Λapm,n,0,0| of the 2D periodic square patch shown in Figure 6.AMC unit cell consisting of a square PEC patch on a dielectric substrate with relative permittivity 3.5. (a) Top view. (b) Side view.
(a)(b)The modes that can be ignored are identified as follows. The right-hand-side vector[Hyi,Hxi]′ is generated for all 2×31×31 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.


## 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 x- and y-directions, is displayed in Figure 1. We denote each homogeneous layer by Ri, i=1,…,n, where n is the number of layers. R0, Rn+1 are the upper and lower half spaces, respectively. The interface between the Ri-1 and Ri regions, denoted by Si, is located at z=zi. 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 Ri has a thickness of hi=zi+1-zi. 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 atz=zn+1;Type B: has a PEC ground plane atz=zn+1.Figure 1
Cross-section of a generic multilayer 2D periodic structure.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)
E⃗i=Exix̂+Eyiŷ+Eziẑ: the incident electric field in Region R0.(b)
E⃗r=Exrx̂+Eyrŷ+Ezrẑ: the reflected electric field in Region R0 if the top conducting pattern at z=z1 were a complete PEC without any apertures. Note that, at z=z1, Exi+Exr=0 and Eyi+Eyr=0.(c)
Mx(Ri,Si): the x-component of the equivalent magnetic currents at the top interface Si of region Ri.(d)
Mx(Ri,Si+1): the x-component of the equivalent magnetic currents at the bottom interface Si+1 of region Ri.(e)
E⃗(Ri)x=Ex(Ri)xx̂+Ey(Ri)xŷ+Ez(Ri)xẑ: the electric field created by Mx(Ri,Si) and Mx(Ri,Si+1) in Region Ri.(f)
E⃗(Ri)s=E⃗(Ri)x+E⃗(Ri)y: the electric field in Region Ri due to all equivalent magnetic currents.(g)
E⃗(Ri)tot=Ex(Ri)totx̂+Ey(Ri)totŷ+Ez(Ri)totẑ: the total electric field in Region Ri.(h)
𝒜(Ri)x=𝒜(Ri,Si)x+𝒜(Ri,Si+1)x: 𝒜(Ri)x is the x-component of the total magnetic vector potential in Region Ri. 𝒜(Ri,Si)x and 𝒜(Ri,Si+1)x are the x-components of the magnetic vector potential created by the top and bottom equivalent magnetic currents, respectively, in Region Ri.(i)
ϵ(Ri), μ(Ri): the permittivity and permeability of the dielectric in region Ri.The total field in regionR0 is given by(1)E⃗(Ri)tot=E⃗(Ri)s+E⃗i+E⃗r,
where the scattered field E⃗(Ri)s is totally due to the magnetic currents at the bottom interface, since there is no top interface for region R0. In other regions,(2)E⃗(Ri)tot=E⃗(Ri)s.
In Type B structures, E⃗(Ri+1)tot=0.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 RegionRi, (1≦i≦n), for example, as shown in Figure 2. The range of Ri is zi+1≦z≦zi. In Figure 2, without any loss of generality, we have assumed that zi+1=0 and zi=h to simplify the explanation. The apertures in the top and bottom interfaces (Si and Si+1) of Region Ri are replaced by equivalent magnetic currents on PECs. The top equivalent magnetic currents (Mx(Ri,Si),My(Ri,Si))  now appear in the aperture areas of Surface Si, just below the PEC at z=zi, while the bottom equivalent magnetic currents (Mx(Ri,Si+1), My(Ri,Si+1)) exist in the apertures of Surface Si+1, just above the bottom PEC at z=zi+1.Figure 2
Typical regionRi. It is assumed that zi+1=0, zi=h.


## 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 Figure4(a), is separated into two regions, Region R0 and Region R1, as shown in Figure 4(b). The metal pattern at z=h 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 Figure4(b), the range of Region R0 is z∈[h,+∞). Its top is open and its bottom has equivalent magnetic currents, Mx(R0) and My(R0), just above a PEC. In Type A structures, the range of Region R1 is z∈(-∞,h] with an open boundary at the bottom. In Type B structures, the range of Region R1 is z∈[0,h] 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 z=h. There exist equivalent magnetic currents, Mx(R1), My(R1), just below this PEC.Next we apply image theory to RegionR0 and remove the PEC surface at z=h, by doubling the magnetic currents. This transforms Region R0 to the “waveguide” shown in Figure 3(a), which has Bloch-Floquet periodic boundaries on four sides, opened at two ends and excited by 2M⃗(R0)=2Mx(R0)x̂+2M⃗y(R0)ŷ at z=h. When image theory is applied to Region R1 of a Type A structure, it is transformed to a similar waveguide, but excited by 2M⃗(R1)=2Mx(R1)x̂+2M⃗y(R1)ŷ at z=h. Region R1 of a Type B structure is transformed to the waveguide shown in Figure 3(c), which is excited by 2M⃗(R1) at z=h and their images at z=(2p+1)h, p=…,-2,-1,1,2,…. 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 z coordinates.Equivalent “waveguides” with Bloch periodic boundaries at four sides, opened at two ends; (a) excited by equivalent magnetic currents atz=h; (b) excited by equivalent magnetic currents at z=0; (c) excited by equivalent magnetic currents and their images at  z=(2p+1)h, p=…,-1,0,1,…; (d) excited by equivalent magnetic currents and their images at z=(2p+1)h and z=(2p)h, where p=…,-1,0,1,….
(a)(b)(c)(d)Figure 4
(a) A unit cell. (b) Separation of a Type B unit cell into two regions. RegionR0 is open at the top and has a PEC at the bottom. Equivalent magnetic currents (shown in light color at z=h) exist just above the bottom PEC. Region R1 has a PEC surface at the top and equivalent magnetic currents (shown in light color at z=h) exist just below the top PEC. For Type A structures, the bottom of Region R1 is open (not shown here). For Type B structures, it is a PEC (as shown).To find the electromagnetic fields produced by magnetic currents in RegionR0 and Region R1 [10], we solve the following wave equations for the magnetic vector potential (𝒜) in the two corresponding waveguides:(3)∇2Ax+k2Ax=-2Mxjωμδ(z-zp+1),(4)∇2Ay+k2Ay=-2Myjωμδ(z-zp+1),(5)E⃗s=-jωμ∇×A⃗,(6)H⃗s=∇×∇×A⃗,
where 𝒜⃗=𝒜xx̂+𝒜yŷ and zp+1=(2p+1)h represent the z coordinate of each layer of equivalent magnetic currents (i.e., excitations in Figures 3(a) and 3(c)). δ(z-zp+1) is infinite at z=zp+1, otherwise it is zero.The general solutions of (3) and (4), the inhomogeneous, second-order partial differential equations, are:(7)A(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)xFX(m)FY(n)FZ(R0,R1)(m,n),(8)A(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)yFX(m)FY(n)FZ(R0,R1)(m,n).
Here, we select FX(m)FY(n) as eigenmodes of the waveguide bounded by Bloch-Floquet periodic boundary conditions. The sum of these modal functions, weighted with coefficients amn(R0,R1)x,y, gives the full-wave solution including all propagating and evanescent waves. The subscript mn represents the sequence of modal functions given below. These modal functions are the same in both regions:(9)FX(m)=e-jkmx,km=2mπa+kx,FY(n)=e-jkny,kn=2nπb+ky,
where kx, ky are the x and y components of the incident wave vector k⃗(R0); a is the lattice constant along the x direction;  b is the lattice constant along the y direction; m and n are integers from -∞ to +∞. In the special case of normal incidence, kx=ky=0.By substituting in the wave equation, one can prove that theFZ(R0)(m,n) in (7) is given by(10)FZ(R0)(m,n)=e-γmn(R0)(z-h),
where(11)(γmn(R0,R1))2=km2+kn2-(k(R0,R1))2,(k(R0,R1))2=ω2ϵ0ϵr(R0,R1)μ0μr(R0,R1).
For a Type A structure, FZ(R1)(m,n) is given by(12)FZ(R1)(m,n)=eγmn(R1)(z-h).
For a Type B structure, FZ(R1)(m,n) is given by(13)FZ(R1)(m,n)=∑p=0+∞[eγmn(R1)(z-(2p+1)h)+e-γmn(R1)(z+(2p+1)h)].
To shorten the mathematical expressions, let us define two new functions, F(R0,R1) and Fd(R0,R1), to represent the 3D field modes in (7) and their derivatives:(14)F(R0,R1)(m,n)=FX(m)FY(n)FZ(R0,R1)(m,n),Fd(R0,R1)(m,n)=FX(m)FY(n)FZd(R0,R1)(m,n),
where(15)FZd(R0,R1)(m,n)=-∂FZ(R0,R1)(m,n)γmn(R0,R1)∂z.
At z=h,(16)FZd(R0)=1,FZd(R1)=-1.
The electric and magnetic fields produced by the magnetic currents in each region are obtained by substituting (7) and (8) into (5) and (6):(17)Ex(R0,R1)x=0,Ey(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)x(jωμ(R0,R1)γmn(R0,R1))Fd(R0,R1)(m,n),Hx(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)x((k(R0,R1))2-(km)2)×F(R0,R1)(m,n),Hy(R0,R1)x=∑m,n=-∞+∞amn(R0,R1)x(-kmkn)F(R0,R1)(m,n),Ex(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)y(-jωμ(R0,R1)γmn(R0,R1))×Fd(R0,R1)(m,n),Ey(R0,R1)y=0,Hx(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)y(-kmkn)F(R0,R1)(m,n),Hy(R0,R1)y=∑m,n=-∞+∞amn(R0,R1)y((k(R0,R1))2-(kn)2)×F(R0,R1)(m,n).


## 2.3. Coupled-Field Expansion for a Unit Cell of a Multilayer Periodic Structure

RegionRi, shown in Figure 2, is the generalized case of Region R0 and Region R1 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 Ri. Its top equivalent magnetic currents Mx(Ri,Si), My(Ri,Si) are doubled when the top PEC at z=zi is removed. The same applies to its bottom equivalent magnetic currents Mx(Ri,Si+1), My(Ri,Si+1) when the bottom PEC at z=zi+1 is removed. In addition, multiple images of these currents appear in an infinite number of surfaces. This transforms Region Ri to a waveguide with Bloch-Floquet periodic boundaries at the four sides, excited by 2Mx(Ri,Si) and 2My(Ri,Si) at z=zi and a series of their images at z=zi+(2p)hi, as well as 2Mx(Ri,Si+1), 2My(Ri,Si+1) at z=zi+1 and a series of their images at z=zi+1+(2p)hi, where p=…,-2,-1,1,2,…. This excitation scheme for a typical Region i is shown in Figure 3(d).As the top and bottom regions are different from a typical region, shown in Figure2, we need to consider four types of “waveguides”, as follows.(1)
Fori=0, Region R0 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, fori=n+1, Region Rn+1 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 Sn+1 is a PEC.(3)
In Type B structures, fori=n, Region Rn leads to the “waveguide” in Figure 3(c) because Surface Sn+1 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 Figure3(d).Note that the last type of “waveguide” is required only for multilayer structures. We exploit the fact that the excitations in Figure3(d) can be realized by shifting and superposing the excitations in Figures 3(a) and 3(b) at an infinite periodic series of z-coordinates.To find the electromagnetic fields in RegionRi produced by magnetic currents on Surfaces Si and Si+1, we solve the following wave equations for the magnetic vector potential (𝒜(Ri)) in the corresponding waveguides:(18)∇2A(Ri)x+(k(Ri))2A(Ri)x=-2Mx(Ri,Si)δ(z-zi-zp)+2Mx(Ri,Si+1)δ(z-zi+1-zp)jωμ(Ri),(19)∇2A(Ri)y+(k(Ri))2A(Ri)y=-2My(Ri,Si)δ(z-zi-zp)+2My(Ri,Si+1)δ(z-zi+1-zp)jωμ(Ri),(20)E⃗(Ri)s=-jωμ(Ri)∇×A⃗(Ri),(21)H⃗(Ri)s=∇×∇×A⃗(Ri),
where 𝒜⃗(Ri)=𝒜(Ri)xx̂+𝒜(Ri)yŷ and zp=(2p)hi, p=…,-2,-1,0,1,2,…. δ(z-zi-(2p)hi) is infinite at z=zi+(2p)hi, otherwise, it is zero. δ(z-zi+1-(2p)hi) is infinite at z=zi+1+(2p)hi, otherwise it is zero. The series of z=zi+zp and z=zi+1+zp   represents the z 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 Figure3(d) can be expressed as a linear superposition of the potentials created by the equivalent magnetic currents at the two interfaces (z=zi and z=zi+1), and their corresponding series of images. So, let us separate them as(22)A(Ri)y=A(Ri,Si+1)y+A(Ri,Si)y,A(Ri)x=A(Ri,Si+1)x+A(Ri,Si)x.
Now, we expand these potentials as(23)A(Ri,Si+1)y=∑m,n=-∞+∞amn(Ri,Si+1)yFX(m)FY(n)FZ(Ri,Si+1)(m,n),(24)A(Ri,Si)y=∑m,n=-∞+∞amn(Ri,Si)yFX(m)FY(n)FZ(Ri,Si)(m,n),(25)A(Ri,Si+1)x=∑m,n=-∞+∞amn(Ri,Si+1)xFX(m)FY(n)FZ(Ri,Si+1)(m,n),(26)A(Ri,Si)x=∑m,n=-∞+∞amn(Ri,Si)xFX(m)FY(n)FZ(Ri,Si)(m,n),
where FX(m) and FY(n) are defined in (9).In (23)–(26), FZ(Ri,Si+1)(m,n) and  FZ(Ri,Si)(m,n) are given by(27)FZ(Ri,Si+1)(m,n)=∑p=0+∞[eγmn(Ri)(z-zi+1-(2p)hi)+e-γmn(Ri)(z-zi+1+(2p)hi)],FZ(Ri,Si)(m,n)=∑p=0+∞[eγmn(Ri)(z-zi-(2p)hi)+e-γmn(Ri)(z-zi+(2p)hi)],
where(28)(γm,n(Ri))2=(km)2+(kn)2-(k(Ri))2,(k(Ri))2=ω2μ(Ri)ϵ(Ri).
To shorten the mathematical expressions, let us define four new functions, F(Ri,Si+1), F(Ri,Si), Fd(Ri,Si+1),  and Fd(Ri,Si)  to represent the 3D field modes in (23) and (26), and their derivatives:(29)F(Ri,Si+1)(m,n)=FX(m)FY(n)FZ(Ri,Si+1)(m,n),F(Ri,Si)(m,n)=FX(m)FY(n)FZ(Ri,Si)(m,n),Fd(Ri,Si+1)(m,n)=FX(m)FY(n)FZd(Ri,Si+1)(m,n),Fd(Ri,Si)(m,n)=FX(m)FY(n)FZd(Ri,Si)(m,n),
where(30)FZd(Ri,Si+1)(m,n)=-1γmn(Ri)⋅∂FZ(Ri,Si+1)(m,n)∂z,FZd(Ri,Si)(m,n)=1γmn(Ri)⋅∂FZ(Ri,Si)(m,n)∂z.
It can be shown that, at z=zi,(31)FZd(Ri,Si+1)=0,FZd(Ri,Si)=1.
Similarly, at z=zi+1,(32)FZd(Ri,Si+1)=1,FZd(Ri,Si)=0.For further compactness, let us define three operators:(33)Γa(Si)m,n[f(x)]=∬apf(x)FX*(m)FY*(n)dxdy,(34)Γp(Si)m,n[f(x)]=∬patchf(x)FX*(m)FY*(n)dxdy,(35)Γu(Si)m,n[f(x)]=∬unitf(x)FX*(m)FY*(n)dxdy,
where f(x) is an arbitrary function. ∬ap, ∬patch, and ∬unit  represent the integration over the aperture areas, metal patch areas, and the whole unit cell area at the interface Si, respectively. The operator Γu(Si)m,n satisfies the following orthogonality:(36)Γu(Si)m,n[FX(m′)FY(n′)]={ab,m′=m,n′=n,0,otherwise.
Because the integration area of a unit cell is equal to the sum of its aperture area and patch area,(37)Γu(Si)m,n[f(x)]=Γa(Si)m,n[f(x)]+Γp(Si)m,n[f(x)].
In the following subsections, various interface continuities are enforced, in order to find the unknown coefficients amn(Ri,Si+1)x,y and amn(Ri,Si)x,y in (23)–(26).Since RegionRi is the generalized case of Region R0 and R1 in Section 2.2, we will enforce the boundary conditions on Surface Si.


## 2.4. Enforcing the Tangential Electric Field to Be Continuous across Apertures in SurfaceSi

Let us consider RegionRi-1, and apply operator (35) to both sides of (48) at z=zi. The result is(38)Γu(Si)m,n[Ex(Ri-1)y]=Γu(Si)m,n[∑m′,n′=-∞+∞(jωμ(Ri-1)γm′n′(Ri-1))⋅(am′n′(Ri-1,Si-1)yFd(Ri-1,Si-1)-am′n′(Ri-1,Si)yFd(Ri-1,Si))∑m′,n′=-∞+∞].
Due to the orthogonality in (36), only the m′=m, n′=n term is nonzero on the right-hand side. Then, considering that Fd(Ri-1,Si-1)|z=zi≡0 and Fd(Ri-1,Si)|z=zi≡1 as given in (32), we have(39)Γu(Si)m,n[Ex(Ri-1)y]=amn(Ri-1,Si)y(-jωμ(Ri-1)γmn(Ri-1))ab.
Therefore,(40)amn(Ri-1,Si)y=Γu(Si)m,n[Ex(Ri-1)y]-jωμ(Ri-1)γmn(Ri-1)ab.
Similarly, we can derive the following expressions by applying the operator (35) to (48) for Region Ri and applying (35) to (53) for Region Ri-1 and Region Ri:
(41)amn(Ri,Si)y=Γu(Si)m,n[Ex(Ri)y]jωμ(Ri)γmn(Ri)ab,(42)amn(Ri-1,Si)x=Γu(Si)m,n[Ey(Ri-1)x]jωμ(Ri-1)γmn(Ri-1)ab,(43)amn(Ri,Si)x=Γu(Si)m,n[Ey(Ri)x]-jωμ(Ri)γmn(Ri)ab.
Let us consider the interface Si in a typical region Ri-1 (i-1≠0,n). The x-component of the total electric field on this plane is given by Ex(Ri-1)tot=Ex(Ri-1)x+Ex(Ri-1)y. Then, as x-directed magnetic currents do not produce an x-directed electric field on the same plane (see (26)), Ex(Ri-1)x=0, and, therefore, Ex(Ri-1)tot=Ex(Ri-1)y at this plane at z=zi. Similarly, we have Ex(Ri)tot=Ex(Ri)y for Region Ri at the same plane at z=zi. 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(44)Γu(Si)m,n[Ex(Ri-1)y]=Γu(Si)m,n[Ex(Ri)y],(45)Γu(Si)m,n[Ey(Ri-1)x]=Γu(Si)m,n[Ey(Ri)x].
In region R0, the x-component of the total electric field is the sum of the x-components of the incident field (Exi), the reflected field (Exr), and the scattered field. However, at the S1 interface, Exi and Exr cancel out due to the PEC, and hence the above relationships are valid also for region R0 at interface S1.Let us substitute (44) into (40)-(41) and make a simple rearrangement which leads to the following relationship between the amn coefficients for the two coupled regions Ri and Ri+1:(46)bmn(Si)y=amn(Ri-1,Si)y(-jωμ(Ri-1)γmn(Ri-1))=amn(Ri,Si)y(jωμ(Ri)γmn(Ri)),
where bmn(Si)y is a new unknown coefficient, which will replace both amn(Ri) and amn(Ri-1) eventually.Similarly, consideringEy in (45) and (42)-(43), we have(47)bmn(Si)x=amn(Ri-1,Si)x(jωμ(Ri-1)γmn(Ri-1))=amn(Ri,Si)x(-jωμ(Ri)γmn(Ri)).
Note that the above two equations establish a simple relationship between the coefficients in Region Ri-1 and those in Region Ri at z=zi. We have exploited this relationship to introduce a new set of unknown coefficients bmn(Si)y and bmn(Si)x, which are common to both regions. Note that the coefficients bmn(Si)x and bmn(Si)y are associated with the ith interface whereas the former coefficients amn(Ri,Sj) 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):(48)Ex(Ri)y=∑m′,n′=-∞+∞bm′n′(Si+1)yFd(Ri,Si+1)(m′,n′)+∑m′,n′=-∞+∞bm′n′(Si)yFd(Ri,Si)(m′,n′),(49)Ey(Ri)y=0,(50)Hx(Ri)y=∑m′,n′=-∞+∞bm′n′(Si+1)y(-km′kn′)F(Ri,Si+1)(m′,n′)-jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)y(-km′kn′)F(Ri,Si)(m′,n′)jωμ(Ri)γm′n′(Ri),(51)Hy(Ri)y=∑m′,n′=-∞+∞bm′n′(Si+1)y((k(Ri))2-kn′2)F(Ri,Si+1)(m′,n′)-jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)y((k(Ri))2-kn′2)F(Ri,Si)(m′,n′)jωμ(Ri)γm′n′(Ri),(52)Ex(Ri)x=0,(53)Ey(Ri)x=∑m′,n′=-∞+∞bm′n′(Si+1)xFd(Ri,Si+1)(m′,n′)+∑m′,n′=-∞+∞bm′n′(Si)xFd(Ri,Si)(m′,n′),(54)Hx(Ri)x=∑m′,n′=-∞+∞bm′n′(Si+1)x((k(Ri))2-km′2)F(Ri,Si+1)(m′,n′)jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)x((k(Ri))2-km′2)F(Ri,Si)(m′,n′)-jωμ(Ri)γm′n′(Ri),(55)Hy(Ri)x=∑m′,n′=-∞+∞bm′n′(Si+1)x(-km′kn′)F(Ri,Si+1)(m′,n′)jωμ(Ri)γm′n′(Ri)+∑m′,n′=-∞+∞bm′n′(Si)x(-km′kn′)F(Ri,Si)(m′,n′)-jωμ(Ri)γm′n′(Ri).


## 2.5. Enforcing the Tangential Electric Field to Be Zero on the Metal Patch Areas in SurfaceSi

As the tangential electric field is zero on the metal patch areas of the interface,Γp(Si)m,n[Ex(Ri-1)y]=0 and Γu(Si)m,n[Ex(Ri-1)y]=Γa(Si)m,n[Ex(Ri-1)y]. Hence, we can replace Γu(Si)m,n in (40) by Γa(Si)m,n as we only consider Ex(Ri-1)y at z=zi. In other words, the integration area in (40) can be reduced to the aperture area:(56)amn(Ri-1,Si)y=Γa(Si)m,n[Ex(Ri-1,Si)y]-jωμ(Ri-1)γmn(Ri-1)ab.
Similarly, one may prove that(57)amn(Ri,Si)y=Γa(Si)m,n[Ex(Ri)y]jωμ(Ri)γmn(Ri)ab,(58)amn(Ri-1,Si)x=Γa(Si)m,n[Ey(Ri-1)x]jωμ(Ri-1)γmn(Ri-1)ab,(59)amn(Ri,Si)x=Γa(Si)m,n[Ey(Ri)x]-jωμ(Ri)γmn(Ri)ab.
Let us  replace amn(Ri-1,Si)y(-jωμ(Ri-1)γmn(Ri-1)) in (56) using bmn(Si)y and expand Ex(Ri-1,Si)y in (56) using (48). We obtain(60)bmn(Si)y=∑m′,n′=-∞+∞bm′n′(Si)yΛap(Si)(m,n,m′,n′)ab.
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(61)bmn(Si)x=∑m′,n′=-∞+∞bm′n′(Si)xΛap(Si)(m,n,m′,n′)ab,
where(62)Λap(Si)(m,n,m′,n′)=Γa(Si)m,n[FX(m′)FY(n′)].
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 SurfaceSi

The totalHy in region R0 is(63)Hy(R0)tot=Hy(R0)x+Hy(R0)y+Hyi+Hyr,
where Hyi and Hyr are the y-components of the magnetic field of incident and reflected waves in region R0. At the interface S1, Hyr=Hyi due to the PEC, and hence(64)Hy(R0)tot=Hy(R0)x+Hy(R0)y+2Hyi.
The total Hy in a typical region, Ri, 1<i<n+1 is(65)Hy(Ri)tot=Hy(Ri)x+Hy(Ri)y.
Now let us apply operator Γa(S1)m,n to (64):(66)Γa(S1)m,n[Hy(R0)tot]=Γa(S1)m,n[Hy(R0)x+Hy(R0)y+2Hyi].
Next, we consider the interface S1 and Region R1 below it. So, we set i=1 in (65) and apply the same operator Γa(S1)m,n:(67)Γa(S1)m,n[Hy(R1)tot]=Γa(S1)m,n[Hy(R1)x+Hy(R1)y].
Since the magnetic field is continuous across the aperture at interface S1, Γa(S1)m,n[Hy(R0)tot]=Γa(S1)m,n[Hy(R1)tot]. So, we subtract (67) from (66). After simply rearranging the sequence of terms, we obtain(68)Γa(S1)m,n[-2Hyi]=Γa(S1)m,n[Hy(R0)x+Hy(R0)y-Hy(R1)x-Hy(R1)y]=(Γu(S1)m,n-Γp(S1)m,n)[Hy(R0)x+Hy(R0)y-Hy(R1)x-Hy(R1)y].
Note that Γa(S1)m,n=Γu(S1)m,n-Γp(S1)m,n due to (37). Similarly the continuity of the tangential component of the magnetic field across the apertures at the other typical interfaces (Si, i>1) leads to(69)0=Γa(Si)m,n[Hy(Ri-1)x+Hy(Ri-1)y-Hy(Ri)x-Hy(Ri)y]=(Γu(Si)m,n-Γp(Si)m,n)×[Hy(Ri-1)x+Hy(Ri-1)y-Hy(Ri)x-Hy(Ri)y].Let us  repeat this process withHx. At interface S1, the relationship is(70)Γa(S1)m,n[-2Hxi]=Γa(S1)m,n[Hx(R0)x+Hx(R0)y-Hx(R1)x-Hx(R1)y]=(Γu(S1)m,n-Γp(S1)m,n)[Hx(R0)x+Hx(R0)y-Hx(R1)x-Hx(R1)y].
For any other typical interface Si, i>1,(71)0=Γa(Si)m,n[Hx(Ri-1)x+Hx(Ri-1)y-Hx(Ri)x-Hx(Ri)y]=(Γu(Si)m,n-Γp(Si)m,n)×[Hx(Ri-1)x+Hx(Ri-1)y-Hx(Ri)x-Hx(Ri)y].Next, we express all terms inside theΓu(Si)m,n-Γp(Si)m,n operator in the right-hand side of (68)–(71) in terms of unknowns bmnx,y. Let us start from Γu(Si)m,n[Hy(Ri-1)x] in (69). Replacing i in (55) with i-1 and then applying the operator, we obtain(72)Γu(Si)m,n[Hy(Ri-1)x]=Γu(Si)m,n[∑m′,n′=-∞+∞bm′n′(Si)x(-km′kn′)F(Ri-1,Si)(m′,n′)jωμ(Ri-1)γm′n′(Ri-1)+∑m′,n′=-∞+∞bm′n′(Si-1)x(-km′kn′)F(Ri-1,Si-1)(m′,n′)-jωμ(Ri-1)γm′n′(Ri-1)].Due to the orthogonality in (36), only the m′=m, n′=n term is nonzero on the right-hand side. Hence,(73)Γu(Si)m,n[Hy(Ri-1)x]=(-kmkn)bmn(Si)xζ(Ri-1,Si)m,nabjωμ(Ri-1)γmn(Ri-1)+(-kmkn)bmn(Si-1)xζ(Ri-1,Si-1)m,nab-jωμ(Ri-1)γmn(Ri-1),
where(74)ζ(Ri-1,Si)m,n=FZ(Ri-1,Si)(m,n)|z=zi,ζ(Ri,Si)m,n=FZ(Ri,Si)(m,n)|z=zi.Nextbmn(Si)x and bmn(Si-1)x are expanded using (61). Similarly, we apply operator Γp(Si)m,n to Hy(Ri)x in Region Ri-1. The result is(75)Γp(Si)m,n[Hy(Ri-1)x]=∑m′,n′=-∞+∞(-km′kn′)bm′n′(Si)xζ(Ri-1,Si)m′,n′(Λ-Λap(Si)(m,n,m′,n′))jωμ(Ri-1)γm′n′(Ri-1)+∑m′,n′=-∞+∞(-km′kn′)bm′n′(Si-1)xζ(Ri-1,Si-1)m′,n′(Λ-Λap(Si)(m,n,m′,n′))-jωμ(Ri-1)γm′n′(Ri-1),
where Λ=ab when m=m′, n=n′, otherwise Λ=0.We similarly expand the other terms in (68)–(71) and express the right-hand side of all these equations in terms of bmnx,y.


## 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 InterfacesSi

The matrix relationship developed in this subsection applies to all intermediate, typical surfacesSi. For Type A structures, they range from i=2,…,n and, for Type B structures, i=2,…,n-1.Let us  define three column vectors to represent unknowns in (69) and (71):(76)B(Si-1)=[B(Si-1)yB(Si-1)x],B(Si)=[B(Si)yB(Si)x],B(Si+1)=[B(Si+1)yB(Si+1)x],
where(77)B(Si-1)y=[…,bm′n′(Si-1)y,…]′;B(Si-1)x=[…,bm′n′(Si-1)x,…]′,B(Si)y=[…,bm′n′(Si)y,…]′;B(Si)x=[…,bm′n′(Si)x,…]′,B(Si+1)y=[…,bm′n′(Si+1)y,…]′;B(Si+1)x=[…,bm′n′(Si+1)x,…]′,
and m′, n′ are from -∞→∞.By substituting expansions (73)-(75), and so forth in (69) and (71), we obtain the linear subsystem for each intermediate layer:(78)0=[ZSi(Si-1)ZSi(Si)ZSi(Si+1)]⋅[B(Si-1)B(Si)B(Si+1)],
where(79)ZSi(Si-1)=[OSi(Si-1)yPSi(Si-1)xQSi(Si-1)yLSi(Si-1)x],ZSi(Si)=[OSi(Si)yPSi(Si)xQSi(Si)yLSi(Si)x],ZSi(Si+1)=[OSi(Si+1)yPSi(Si+1)xQSi(Si+1)yLSi(Si+1)x].
The above expression describes the relationship between the unknowns for three adjacent surfaces. In this expression, OSi(Si-1)y, OSi(Si)y, OSi(Si+1)y, PSi(Si-1)x, PSi(Si)x, PSi(Si+1)x, QSi(Si-1)y, QSi(Si)y, QSi(Si+1)y, LSi(Si-1)x, LSi(Si)x, and LSi(Si+1)x are twelve submatrices. Their elements are(80)OSi(Si-1)y=τ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′-τ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,PSi(Si-1)x=-ξ(Ri-1)m,mζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′+ξ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,QSi(Si-1)y=ξ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,LSi(Si-1)x=-φ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′+φ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,OSi(Si)y=(τ(Ri)m,nζ(Ri,Si)m,n-τ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′-(τ(Ri)m′,n′ζ(Ri,Si)m′,n′-τ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,PSi(Si)x=-(ξ(Ri)m,nζ(Ri,Si)m,n-ξ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′+(ξ(Ri)m′,n′ζ(Ri,Si)m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,QSi(Si)y=(ξ(Ri)m,nζ(Ri,Si)m,n-ξ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′-(ξ(Ri)m′,n′ζ(Ri,Si)m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,LSi(Si)x=-(φ(Ri)m,nζ(Ri,Si)m,n-φ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′+(φ(Ri)m′,n′ζ(Ri,Si)m′,n′-φ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,OSi(Si+1)y=-τ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′+τ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,PSi(Si+1)x=ξ(Ri)m,mζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′-ξ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,QSi(Si+1)y=-ξ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′+ξ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,LSi(Si+1)x=φ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′-φ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,
where Λp(Si)m,n,m′,n′=Λ-Λap(Si)m,n,m′,n′, and ξ(Ri)m,n, τ(Ri)m,n, φ(Ri)m,n are given by(81)ξ(Ri)m,n=(-kmkn)jωμ(Ri)γmn(Ri),τ(Ri)m,n=((k(Ri))2-kn2)jωμ(Ri)γm,n(Ri),φ(Ri)m,n=((k(Ri))2-km2)jωμ(Ri)γm,n(Ri).


### 2.7.2. Last Interface (Sn+1)

When the derivation in the above subsection is repeated for the last interfaceSn+1 of a Type A structure, that is, Sn+1, the result is(82)0=[ZSn+1(Sn)ZSn+1(Sn+1)]⋅[B(Sn)B(Sn+1)].
This involves only two unknown vectors because there are no magnetic currents below this surface.By manipulating the above equation,B(Sn+1) can be expressed in terms of B(Sn):(83)B(Sn+1)=[ZSn+1(Sn+1)]-1[-ZSn+1(Sn)B(Sn)].
Let us now define(84)Φn+1=-[ZSn+1(Sn+1)]-1ZSn+1(Sn),
then (83) becomes(85)B(Sn+1)=Φn+1B(Sn).
Note that now we have a relationship connecting unknowns for two adjacent interfaces.In a TypeB structure, the last surface with magnetic currents is Sn. The above derivation is valid for this case too. In this case, B(Sn+1)=0; the relationship between B(Sn) and B(Sn-1) is similar to (85).


### 2.7.3. Second-Last InterfaceSn

As an example, let us now consider the next interface,Sn. The objective here is to derive a relationship between B(Sn) and B(Sn-1), by exploiting (78) and (85). Let us first substitute i=n in (78) and rewrite it as(86)ZSn(Sn)B(Sn)+ZSn(Sn+1)B(Sn+1)=-ZSn(Sn-1)B(Sn-1).
Now, we substitute (85) to replace B(Sn+1). The result is(87)[ZSn(Sn)+ZSn(Sn+1)Φn+1]B(Sn)=-ZSn(Sn-1)B(Sn-1).
After rearrangement, we get the following relationship between B(Sn) and B(Sn-1):(88)B(Sn)=ΦnB(Sn-1),
where(89)Φn=-[ZSn(Sn)+ZSn(Sn+1)Φn+1]-1⋅ZSn(Sn-1).
In a Type B structure, the second-last surface with magnetic currents is Sn-1. The above derivation is valid for this case too.


### 2.7.4. Intermediate InterfaceSi

Now, we revisit a typical intermediate surfaceSi. As shown above for the i=n example, it is possible to express the unknown coefficients B(Si) (i=2,…,n+1) in terms of B(Si-1) by repeating the same steps. The result is(90)B(Si)=ΦiB(Si-1),
where(91)Φi=-[ZSi(Si)+ZSi(Si+1)Φi+1]-1⋅ZSi(Si-1).


### 2.7.5. Top InterfaceS1

This interface is different due to the existence of the incident field in RegionR0. After expanding (68) and (70), we obtain(92)H(S1)=[ZS1(S1)ZS1(S2)]⋅[B(S1)B(S2)],
where H(S1) is defined as(93)H(S1)=[Hy(S1)iHx(S1)i],Hy(S1)i=[…,Γam,n[-2Hyi],…]′,Hx(S1)i=[…,Γam,n[-2Hxi],…]′.
Using (90), B(S2) in (92) can be replaced by B(S1). The result is(94)H(S1)=[ZS1(S1)+ZS1(S2)Φ2]⋅[B(S1)].
Note that the only unknown here is B(S1), which can be obtained by solving this matrix equation. Once B(S1) is found, all other unknowns can be found, if necessary, step by step, following (90).When calculating the reflection characteristics of multilayer structures, onlyB(S1) is required. In fact, if all modes other than the m=0, n=0 are evanescent, one needs only b00(S1) to determine the surface reflection coefficients.


## 2.7.1. Intermediate InterfacesSi

The matrix relationship developed in this subsection applies to all intermediate, typical surfacesSi. For Type A structures, they range from i=2,…,n and, for Type B structures, i=2,…,n-1.Let us  define three column vectors to represent unknowns in (69) and (71):(76)B(Si-1)=[B(Si-1)yB(Si-1)x],B(Si)=[B(Si)yB(Si)x],B(Si+1)=[B(Si+1)yB(Si+1)x],
where(77)B(Si-1)y=[…,bm′n′(Si-1)y,…]′;B(Si-1)x=[…,bm′n′(Si-1)x,…]′,B(Si)y=[…,bm′n′(Si)y,…]′;B(Si)x=[…,bm′n′(Si)x,…]′,B(Si+1)y=[…,bm′n′(Si+1)y,…]′;B(Si+1)x=[…,bm′n′(Si+1)x,…]′,
and m′, n′ are from -∞→∞.By substituting expansions (73)-(75), and so forth in (69) and (71), we obtain the linear subsystem for each intermediate layer:(78)0=[ZSi(Si-1)ZSi(Si)ZSi(Si+1)]⋅[B(Si-1)B(Si)B(Si+1)],
where(79)ZSi(Si-1)=[OSi(Si-1)yPSi(Si-1)xQSi(Si-1)yLSi(Si-1)x],ZSi(Si)=[OSi(Si)yPSi(Si)xQSi(Si)yLSi(Si)x],ZSi(Si+1)=[OSi(Si+1)yPSi(Si+1)xQSi(Si+1)yLSi(Si+1)x].
The above expression describes the relationship between the unknowns for three adjacent surfaces. In this expression, OSi(Si-1)y, OSi(Si)y, OSi(Si+1)y, PSi(Si-1)x, PSi(Si)x, PSi(Si+1)x, QSi(Si-1)y, QSi(Si)y, QSi(Si+1)y, LSi(Si-1)x, LSi(Si)x, and LSi(Si+1)x are twelve submatrices. Their elements are(80)OSi(Si-1)y=τ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′-τ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,PSi(Si-1)x=-ξ(Ri-1)m,mζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′+ξ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,QSi(Si-1)y=ξ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,LSi(Si-1)x=-φ(Ri-1)m,nζ(Ri-1,Si-1)m,nΛap(Si)m,n,m′,n′+φ(Ri-1)m′,n′ζ(Ri-1,Si-1)m′,n′Λp(Si)m,n,m′,n′,OSi(Si)y=(τ(Ri)m,nζ(Ri,Si)m,n-τ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′-(τ(Ri)m′,n′ζ(Ri,Si)m′,n′-τ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,PSi(Si)x=-(ξ(Ri)m,nζ(Ri,Si)m,n-ξ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′+(ξ(Ri)m′,n′ζ(Ri,Si)m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,QSi(Si)y=(ξ(Ri)m,nζ(Ri,Si)m,n-ξ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′-(ξ(Ri)m′,n′ζ(Ri,Si)m′,n′-ξ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,LSi(Si)x=-(φ(Ri)m,nζ(Ri,Si)m,n-φ(Ri-1)m,nζ(Ri-1,Si)m,n)Λap(Si)m,n,m′,n′+(φ(Ri)m′,n′ζ(Ri,Si)m′,n′-φ(Ri-1)m′,n′ζ(Ri-1,Si)m′,n′)Λp(Si)m,n,m′,n′,OSi(Si+1)y=-τ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′+τ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,PSi(Si+1)x=ξ(Ri)m,mζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′-ξ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,QSi(Si+1)y=-ξ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′+ξ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,LSi(Si+1)x=φ(Ri)m,nζ(Ri,Si+1)m,nΛap(Si)m,n,m′,n′-φ(Ri)m′,n′ζ(Ri,Si+1)m′,n′Λp(Si)m,n,m′,n′,
where Λp(Si)m,n,m′,n′=Λ-Λap(Si)m,n,m′,n′, and ξ(Ri)m,n, τ(Ri)m,n, φ(Ri)m,n are given by(81)ξ(Ri)m,n=(-kmkn)jωμ(Ri)γmn(Ri),τ(Ri)m,n=((k(Ri))2-kn2)jωμ(Ri)γm,n(Ri),φ(Ri)m,n=((k(Ri))2-km2)jωμ(Ri)γm,n(Ri).


## 2.7.2. Last Interface (Sn+1)

When the derivation in the above subsection is repeated for the last interfaceSn+1 of a Type A structure, that is, Sn+1, the result is(82)0=[ZSn+1(Sn)ZSn+1(Sn+1)]⋅[B(Sn)B(Sn+1)].
This involves only two unknown vectors because there are no magnetic currents below this surface.By manipulating the above equation,B(Sn+1) can be expressed in terms of B(Sn):(83)B(Sn+1)=[ZSn+1(Sn+1)]-1[-ZSn+1(Sn)B(Sn)].
Let us now define(84)Φn+1=-[ZSn+1(Sn+1)]-1ZSn+1(Sn),
then (83) becomes(85)B(Sn+1)=Φn+1B(Sn).
Note that now we have a relationship connecting unknowns for two adjacent interfaces.In a TypeB structure, the last surface with magnetic currents is Sn. The above derivation is valid for this case too. In this case, B(Sn+1)=0; the relationship between B(Sn) and B(Sn-1) is similar to (85).


## 2.7.3. Second-Last InterfaceSn

As an example, let us now consider the next interface,Sn. The objective here is to derive a relationship between B(Sn) and B(Sn-1), by exploiting (78) and (85). Let us first substitute i=n in (78) and rewrite it as(86)ZSn(Sn)B(Sn)+ZSn(Sn+1)B(Sn+1)=-ZSn(Sn-1)B(Sn-1).
Now, we substitute (85) to replace B(Sn+1). The result is(87)[ZSn(Sn)+ZSn(Sn+1)Φn+1]B(Sn)=-ZSn(Sn-1)B(Sn-1).
After rearrangement, we get the following relationship between B(Sn) and B(Sn-1):(88)B(Sn)=ΦnB(Sn-1),
where(89)Φn=-[ZSn(Sn)+ZSn(Sn+1)Φn+1]-1⋅ZSn(Sn-1).
In a Type B structure, the second-last surface with magnetic currents is Sn-1. The above derivation is valid for this case too.


## 2.7.4. Intermediate InterfaceSi

Now, we revisit a typical intermediate surfaceSi. As shown above for the i=n example, it is possible to express the unknown coefficients B(Si) (i=2,…,n+1) in terms of B(Si-1) by repeating the same steps. The result is(90)B(Si)=ΦiB(Si-1),
where(91)Φi=-[ZSi(Si)+ZSi(Si+1)Φi+1]-1⋅ZSi(Si-1).


## 2.7.5. Top InterfaceS1

This interface is different due to the existence of the incident field in RegionR0. After expanding (68) and (70), we obtain(92)H(S1)=[ZS1(S1)ZS1(S2)]⋅[B(S1)B(S2)],
where H(S1) is defined as(93)H(S1)=[Hy(S1)iHx(S1)i],Hy(S1)i=[…,Γam,n[-2Hyi],…]′,Hx(S1)i=[…,Γam,n[-2Hxi],…]′.
Using (90), B(S2) in (92) can be replaced by B(S1). The result is(94)H(S1)=[ZS1(S1)+ZS1(S2)Φ2]⋅[B(S1)].
Note that the only unknown here is B(S1), which can be obtained by solving this matrix equation. Once B(S1) is found, all other unknowns can be found, if necessary, step by step, following (90).When calculating the reflection characteristics of multilayer structures, onlyB(S1) is required. In fact, if all modes other than the m=0, n=0 are evanescent, one needs only b00(S1) 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 m and n as [-15,…,0,…,15], then we have a linear system with 2×31×31≈2000 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 as1/R1-ν in general, where R is the distance to the singular centres, such as aperture corners, edges and tips, and ν is the singularity exponent (e.g., ν=0.5 for a singularity along the edge of an infinite thin metal patch). Its spatial spectrum is given by 2Ψ(ν)(k)-νcos(νπ/2), 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 Ψ(ν)=1.7725 according to [11]. This spectrum is effectively truncated when we truncate the modes by limiting the range of m and n. For this type of singularity in a 10mm unit cell, when the number of modes selected (i.e., maxm,n) is 7, its truncation error is less than 2.4%. There should be at least 15 modes to have an error less than 1%. Hence, the range of m and n is limited to [-15,…,0,…,15] for the calculations presented here. This would normally lead to 31×31 unknowns, or altogether 2×31×31 in both x- and y-directions.However, not all the2×31×31 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 [Hyi,Hxi]′. It can be shown that the ((m,n)th) elements of Hyi and Hxi are proportional to Λapm,n,0,0, which is determined by the geometry of the conducting pattern. The |Λapm,n,0,0| of the periodic square patch, depicted in Figure 6, is shown in Figure 5. We note that |Λapm,n,0,0| along diagonal directions decay much faster than those along the axes, that is, |Λapm,0,0,0| or |Λap0,n,0,0|. It is significant only along the axes and around the origin. This makes the on-axis unknowns bm0 and b0n dominant over off-axis unknowns. Hence, most of the off-axis bmn (i.e., m≠0 and n≠0) can be ignored in the linear system without significantly affecting accuracy.Figure 5
|Λapm,n,0,0| of the 2D periodic square patch shown in Figure 6.AMC unit cell consisting of a square PEC patch on a dielectric substrate with relative permittivity 3.5. (a) Top view. (b) Side view.
(a)(b)The modes that can be ignored are identified as follows. The right-hand-side vector[Hyi,Hxi]′ is generated for all 2×31×31 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 thanm=0 and n=0 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 (RTE and TTE) for the m=0, n=0 mode are obtained as(95)RTE=(b00yux+b00xuy+ErExiux+Eyiuy)atz=h=(b00yux+b00xuyExiux+Eyiuy-1)atz=h,TTE=(b00yux+b00xuyExiux+Eyiuy)atz=h,
where (ux,uy,uz) are the three coordinate components of the unit vector û, which is perpendicular to both the propagation direction vector k̂ and the surface normal direction ẑ:(96)û=k̂×ẑ.
Similarly, for TM-polarized incident waves, the reflection and transmission coefficients are(97)RTM=-(Hysuy+Hxsux+HrHxiux+Hyiuy)atz=h=-(Hysuy+HxsuxHxiux+Hyiuy+1)atz=h,TTM=-η(2)η(1)(Hysuy+HxsuxHxiux+Hyiuy)atz=h,
where(98)Hxs=(-kxky)b00y-jωμ(1)γ00(1)+((k(1))2-(kx)2)b00xjωμ(1)γ00(1),Hys=((k(1))2-(ky)2)b00y-jωμ(1)γ00(1)+(-kxky)b00xjωμ(1)γ00(1),η(1,2)=μ(1,2)ϵ(1,2).
For convenience and without any loss of generality, we chose to calculate the reflection and transmission coefficients at x, y=0. For the special case of normal incidence,(99)R=(Exr+b00yExi)atz=h=(b00yExi-1)atz=h,T=(b00yEi)atz=h.
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 R/T calculation in HFSS cannot be used to calculate RTE,TM and TTE,TM. 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 b00y and b00x by integrating the total field over the aperture area, according to(100)b00y=Γa00[Extot|z=h]ab=∬apExtot|z=hdxdyab,b00x=Γa00[Eytot|z=h]ab=∬apEytot|z=hdxdyab,
and then calculate RTE,TM and TTE,TM 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 compareRTE,TM and TTE,TM 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 (θ=180°, ϕ=90°), has been calculated. As shown in Figure 7, excellent convergence was achieved with 593 Bloch-Floquet modes in the x-direction and the same number of Bloch-Floquet modes in the y-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.Figure 7
CFEM convergence for the square patch AMC in Figure6. Reflection phase angle at normal incidence θ=180° is shown for different CFEM modes and HFSS methods. Substrate thickness is 1.57mm. “tet” represents the number of tetrahedrons.To assess the effect of the substrate thickness on the convergence, we reduced the substrate thickness in Figure6 to 0.5 mm. Figure 8 shows that the method still converges with 593×2 modes. This is because we analytically describe the effect of the substrate thickness in our model.Figure 8
CFEM convergence for the square patch AMC in Figure6. Reflection phase angle at normal incidence (θ=180°) is shown for different CFEM modes and HFSS methods. Substrate thickness is 0.5mm. “tet” represents the number of tetrahedrons.The TE and TM reflection phase for oblique incidence (θ=150°, ϕ=90°) has been obtained with 593×2 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 ±90°, the AMC band for TM incidence is shifted upwards compared to TE incidence.Figure 9
The reflection phase of the square patch AMC in Figure6, for TE oblique incidence at θ=150°, ϕ=90° is shown and compared with the results obtained using HFSS (near-field and far-field methods).Figure 10
The reflection phase of the square patch AMC in Figure6, for TM oblique incidence at θ=150°, ϕ=90° is shown and compared with the results obtained using HFSS (near-field and far-field methods).From Figures7–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 Figure11, 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 a=3mm, b=4mm, w=0.25mm, w1=0.25mm, w2=0.25mm, l=3.75mm, l1=3.75mm, l2=1.5mm; the thickness of the substrate is t=1.6mm; the relative permittivity of the substrate is 4.4. An incident wave linearly polarized at 45° to the x-axis has equal x- and y-components (Exi, Eyi). When all the modes other than m=0, n=0 are evanescent, the reflection coefficients for Ex and Ey polarizations are(101)Rx=(b00(S1)yExi-1)|z=z1,Ry=(b00(S1)xEyi-1)|z=z1.
This circular polarizer was designed such that the reflection coefficient magnitudes (|Rx| and |Ry|) are close to 0dB over a wide band (9.5–14.5 GHz) and the difference between the phases of the reflection coefficients (∠Rx-∠Ry) 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.Geometry of a circular polarizer. (a) Metallic pattern on the top surface. (b) Cross-section. (c) Metallic pattern on the bottom surface.
(a)(b)(c)Figure 12
Reflection phases of two orthogonal field components of the circular polarizer, calculated using CFEM and HFSS.


## 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 (1.57mm) and thin (0.5mm) 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.



---

*Source: 170963-2012-05-10.xml*