Improved Generalized Single-Source Tangential Equivalence Principle Algorithm with Contact-Region Modeling Method for Array Structures

An improved generalized single-source tangential equivalence principle algorithm (GSST-EPA) is proposed for analyzing array structures with connected elements. In order to use the advantages of GSST-EPA, the connected array elements are decomposed and computed by a contact-region modeling (CRM) method, which makes that each element has the same meshes. The unknowns of elements can be transferred onto the equivalence surfaces by GSST-EPA. The scattering matrix in GSST-EPA needs to be solved and stored only once due to the same meshes for each element. The shift invariant of translation matrices is also used to reduce the computation of near-field interaction. Furthermore, the multilevel fast multipole algorithm (MLFMA) is used to accelerate the matrix-vector multiplication in the GSST-EPA. Numerical results are shown to demonstrate the accuracy and efficiency of the proposed method.


Introduction
Surface integral equation (SIE) methods have been widely used in the simulation of complex electromagnetic phenomena in practical engineering.Application of SIE methods always leads to a dense matrix which limits its applicability.Many fast solvers have been developed to reduce the memory requirement and CPU time, such as fast multipole methods [1,2], FFT-based methods [3,4], low-rank-based methods [5][6][7][8], and direct solvers [9,10].However, for many multiscale electromagnetic problems, the dimension of the matrix is large and the matrix is ill-conditioned, which are still the challenges for the fast solvers.It is desirable for more efficient and robust SIE solutions to deal with those challenges.
In recent years, domain decomposition (DD) methods have been attracting much attention and regarded as effective techniques to solve multiscale systems.The DD methods rely on decomposing the problem into several smaller subdomains and solving each subdomain independently.In this manner, the dimension of the matrix can be reduced and the condition number of the matrix can be improved.
Due to these advantages, the DD methods have been combined with finite element (FE) methods [11,12] and finite difference (FD) methods [13,14] successfully.To extend the DD methods to SIE methods, many studies have been proposed.Among these studies, several works are worth mentioning.The overlapping IE-DDM [15] needs to construct buffer regions.The IE-DDM in [16] uses artificial interfaces to do the decomposition.A nonoverlapping IE-DDM [17] has been proposed by using the explicit boundary condition.The discontinuous Galerkin integral equation (IEDG) [18] has been formulated by employing an interior penalty term, which shows great flexibility to geometries.The equivalence principle algorithm (EPA) [19] is a kind of DD method based on Huygens' principle.By using the virtual equivalence surfaces that enclose the subdomains, the EPA can isolate the subdomain from others.Then, the original unknowns of subdomains are transferred to the unknowns on the equivalence surfaces, which is different from the above methods.The linear embedding via Green's operator (LEGO) [20] and generalized transition matrix (GTM) algorithm [21] are similar with EPA.However, the connected array structures have not been dealt with in those methods.
In this paper, the improved generalized single-source tangential equivalence principle algorithm (GSST-EPA) is proposed to address the array structures with connected elements.Although the GSST-EPA has been developed to solve the connected structures [22], it is inefficient for the array structures.For the array structures with connected elements, if each element is decomposed as a subdomain and enclosed by the same equivalence surfaces, then the meshes of elements are different from each other.Therefore, the scattering matrix of each equivalence surface is different and needs to be solved and stored one by one, which leads to huge computational time and memory requirement.To overcome the above problem and make full use of the GSST-EPA, the contact-region modeling (CRM) [23] method is applied to make the decomposition for the array structures.Unlike many other DD methods, the CRM method automatically contains the boundary conditions for the contact surface between subdomains, so the additional transmission conditions are unnecessary.There are two major benefits of the proposed scheme: (a) Each element has the same meshes, so the scattering matrix of each equivalence surface is the same with each other and needs to be solved and stored only once, and (b) the shift invariant of translation matrices can be applied to reduce the computation of coupling between neighbor elements.Finally, the proposed method is accelerated by the multilevel fast multipole algorithm (MLFMA) [2].

The Domain Decomposition Based on the CRM Method
Without loss of generality, consider the scattering problem of a three-dimensional dielectric domain with the media ε 1 , μ 1 , as illustrated in Figure 1(a).The exterior region Ω 0 is the free space.The dielectric structure is illuminated by an incident plane wave E inc , H inc .This problem can be solved by the commonly used PMCHWT formulation [24].In this work, the CRM method is used for the domain partitioning.As shown in Figure 1(b), the dielectric domain Ω is decomposed into two nonoverlapping subdomains such that Ω = Ω 1 ∪ Ω 2 .The surfaces of the subdomains Ω 1 and Ω 2 with outward-directed unit normal n1 and n2 are denoted by S 1 and S 2 , respectively.As the distance δ → 0, the solution of Figure 1(b) is the same as that of Figure 1(a) due to the same physics.Next, the variational PMCHWT formulation based on the CRM method is derived.
By using the extinction theorem, the total fields on the inner boundary of subdomain Ω 1 are equal to zero in the exterior equivalence situation.J 1 and M 1 are the equivalent electric and magnetic currents defined on the outer side of the surface S 1 .The EFIE and MFIE for subdomain Ω 1 can be expressed as where η 0 is the wave impedance in free space.S − 1 denotes the inner side of S 1 .The operators ℒ and K are defined as where X is J or M, i denotes the media domain, and S j and S k are the field surface and source surface, respectively.k i stands for the wave number in domain i.G i is Green's function defined as PV stands for the Cauchy principle value integration.For the residue term ± 1/2 X r × nk , if the observation point r resides on the outer side of S k , then "+" is chosen.Otherwise, "−" is taken if r resides on the inner side of S k , while if r is not on S k , the residue term should be equal to zero.Similarly, the EFIE and MFIE for subdomain Ω 2 in the exterior equivalence situation can be written as The scattering problem of a three-dimensional dielectric object.(b) The object decomposed into two subdomains using the contact-region modeling method.

2
International Journal of Antennas and Propagation For the interior equivalence situation, the total fields are zero on the outer boundary of Ω 1 .So the EFIE and MFIE for the Ω 1 are written as where S + 1 denotes the outer side of S 1 .In the domain Ω 2 , the EFIE and MFIE are Combining (1) and ( 8), the PMCHWT E-field equation for surface S 1 can be derived as The H-field equation can be obtained by combining ( 2) and ( 9) as Similarly, the PMCHWT E-field and H-field equations for surface S 2 can be derived by combining ( 6) and ( 11) and ( 7) and ( 12), as follows: Finally, the matrix form of these four PMCHWT equations is written as where η 0 is used to balance the magnitude of operators and make the matrix better conditioned.

Boundary Conditions at the Contact Surface in CRM
It is worth mentioning that the additional transmission conditions are unnecessary for the contact surface, because the continuity of the currents at the contact surface is contained in the above PMCHWT equations automatically.
To prove this fact, a simple model is considered as shown in Figure 2. The surfaces of domains Ω 1 and Ω 2 are divided into four parts as S 1 = S d1 ∪ S c1 and S 2 = S d2 ∪ S c2 .In the exterior equivalence situation, the MFIE on the contact surfaces S c1 and S c2 are 3 International Journal of Antennas and Propagation where Because the surfaces S c1 and S c2 touch each other, the outer side of S c1 is the same as that of S c2 .Therefore, the terms , and H inc S c1 in (17) are equal to the terms , and H inc S c2 in (18).Since the K operator does not have the residue term, the terms K 0 S c1 S d1 J d1 and K 0 S c1 S d2 J d2 are the same as K 0 S c2 S d1 J d1 and K 0 S c2 S d2 J d2 .By subtracting (18) from (17), the following can be obtained: The definition of the K operator can be found in (4).It can be seen that the Cauchy principle value integration in K 0− S c1 S c1 is the same as that in K 0+ S c2 S c1 , while the residue terms are different.The terms K 0+ S c1 S c2 and K 0− S c2 S c2 have similar relations.Therefore, (20) can be rewritten as Since nc1 = −n c2 , then (21) can be derived as It can be seen that the continuity of the electric currents at the contact surface is included in the MFIE.Similarly, the continuity of the magnetic currents at the contact surface can be derived from the EFIE equations.

GSST-EPA with the CRM Method
Consider the scattering of a dielectric object in free space.As shown in Figure 3, the object is decomposed into N-connected subdomains using the CRM method.Then, each subdomain is enclosed by an equivalence surface.The basic idea of GSST-EPA is to transform the original scattering problem into a new equivalent problem with a scattering operator and a translation operator.Next, the definitions of the scattering and translation operators are illustrated in the following.
The scattering operator can be derived from a dielectric object enclosed by an equivalence surface with three steps as shown in Figure 4.The three steps are outside-in propagation, solving for the currents on the object, and inside-out propagation.The equivalence current on the equivalence surface generates the same scattered field as that generated by original currents on the object.Therefore, the scattering of a dielectric object can be characterized by the scattering operator as Figure 2: The model for the continuity of the currents at the contact surface.where subscripts l and d denote the equivalence surface and dielectric object separately.TEH means one type of the expressions between the incident current J i l and the scattering current J s l , which can be found in [25].The dielectric object is solved by PMCHWT equations.The scattering operator can capture the full physics of the object, so the inside object can be removed.
The coupling between objects via the equivalence surfaces can be described by a translation operator.The translation operator has two different forms: one is suitable for nonadjacent cases, and the other is suitable for adjacent cases.For the two nonadjacent subdomains, the translation operator is quite simple.As shown in Figure 5, the subdomains l and k are enclosed by equivalence surfaces ES l and ES k separately.The coupling from ES k to ES l is the fields produced by the equivalent scattered current J s k on ES k, which can be expressed via a translation operator T lk , l ≠ k, as Therefore, the equivalence scattered current J s k on ES l is given by The expression of the translation operator T kl is similar as T lk given above.
For the two adjacent subdomains, the translation operator can be derived by the source reconstruction method (SRM) [26][27][28].As shown in Figure 6, the inside object has been decomposed into two subdomains with the distance δ → 0. The coupling from ES k to ES l can be derived with several steps.Firstly, suppose the J s k on ES k is already known, then the equivalent electric current J d k and magnetic current M d k on the surface of dielectric k can be obtained by using the SRM as where • † is the pseudoinverse.Secondly, once the equivalent currents on the surface of dielectric k is obtained, the coupling from the dielectric k to the dielectric l can be solved using the CRM method as Thirdly, the equivalent current produced by the incident fields E i d l d k and H i d l d k using the equivalence principle can be written as n l ˆnl ˆnl Figure 4: A dielectric object enclosed by an equivalence surface with three steps.
Finally, by combing ( 27), (28), and ( 29), the equivalence scattered current J s l is equal to where the translation operator T SRM lk is defined as the multiplication of several operators as Since the scattering operator and translation operator have been derived, the GSST-EPA equations for the object in Figure 3 can be written as where

33
ℐ is the identity operator.The matrix form can be obtained by expanding the currents with the Rao-Wilton-Glisson (RWG) basis function and using Galerkin's method to discretize the equations.

Computational Considerations
In this paper, the periodic array structures will be solved by the GSST-EPA.If the array elements are connected with each other, then the CRM method can be used to decompose the elements.Next, each element is enclosed by an equivalence surface.Due to the same geometry and meshes of each element and equivalence surface, the scattering operator S TEH ii is also the same with each other.Therefore, the scattering matrix needs to be solved and stored only once, which can reduce the computational time and memory requirement.
Since the array structures have a lot of elements, the MLFMA is used to accelerate the matrix-vector multiplication in the GSST-EPA.It is notable that the translation operators have two different forms.The T TEH ij operator denotes the coupling between nonadjacent elements which can be calculated by MLFMA, while the T SRM ij operator denotes the near-field coupling between adjacent elements that should be solved directly.Due to characteristic of the periodic array structure, the translation matrices are shift invariant which can be used to reduce the computation of 6 International Journal of Antennas and Propagation near-field coupling.As shown in Figure 7, it is easy to find that the eight kinds of near-interaction mechanism between nine elements is a typical representation of different near interactions between all elements.In fact, the four corner elements in the array only have three kinds of near interactions.The edge elements have five kinds of near interactions.The other elements have eight kinds of near interactions.According to the shift invariant, these eight kinds of T SRM ij can be solved and stored only once, which can be used to reduce the calculation of any adjacent interactions.

Numerical Results
In this section, the performance of the proposed method is studied via several numerical experiments.All of the simulations are performed on PC with Intel Core i7-2760 2.4 GHz CPU and 32.0 GB memory.
Firstly, the accuracy of the CRM is investigated.The scattering of a dielectric sphere with ε r = 2 0 and μ r = 1 0 in free space is considered as shown in Figure 8(a).The sphere is partitioned into two subdomains (two hemispheres) with a fictitious contact surface by CRM as shown in Figure 8(b).The excitation is an x-polarized plane wave propagating into the negative ẑ direction at 300 MHz.The curvilinear Rao-Wilton-Glisson (CRWG) [29] is used to discretize the sphere surface.The bistatic radar cross section (RCS) of HH polarization is given in Figure 9.It can be seen that the of CRM agrees very well with Mie series solution, which also proves that the decomposition does not change the physics of the original problem.
The second example is used to investigate the accuracy and convergence of the GSST-EPA with CRM.As shown in Figure 10, the dimension of the dielectric slab is 1 5 × 1 5 × 0 3 m 3 .The relative permittivity and permeability of the slab are ε r = 2 0 and μ r = 1 0. The excitation is an x-polarized plane wave propagating into the negative ẑ direction at 300 MHz.The slab is decomposed into 3 × 3, 4 × 4, and 5 × 5 subdomains separately by CRM as shown in Figure 11.The subdomains have the same size as each other.Then, each subdomain is enclosed by an equivalence surface with the same size.The sizes of each equivalence surface are 0 7 × 0 7 × 0 5 m 3 , 0 6 × 0 6 × 0 5 m 3 , and 0 5 × 0 5 × 0 5 m 3 for the three different decompositions separately, which means  the distance between the subdomain and the equivalence surface is 0 1λ.Since each subdomain and each equivalence surface have the same size, the same meshes for each subdomain and each equivalence surface can be constructed.Therefore, the scattering matrix in GSST-EPA only needs to be calculated only once.The bistatic RCS of HH polarization is given in Figure 12, which shows good accuracy of GSST-EPA with CRM using different decomposition schemes.Figure 13 shows the comparison of the number among the three different decompositions in GSST-EPA.The generalized minimal residual (GMRES) method with the tolerance 10 −6 is used to solve the final equations.With the increase in subdomains, the iteration number only grows  8 International Journal of Antennas and Propagation from 20 to 25, which shows the convergence stability of the method, while the PMCHWT equations converge to 10 −6 with 1219 steps.
The third example is a 3 × 3 bowtie antenna array located in the x-y plane.The dimension of a bowtie antenna element is shown in Figure 14(a).The yellow part is the PEC patch, while the red part is the dielectric substrate.The thickness of the antenna substrate is 0.03 m with ε r = 2 2 and μ r = 1 0. The elements are connected with each other.The array is divided into 9 subdomains by CRM, which means that each element is a subdomain.Then, each element is by an equivalence surface as shown in Figure 14(b).The excitation is an x-polarized plane wave propagating into the negative ẑ direction at 1.0 GHz.Because the element is a composite target, the EFIE-PMCHWT is used to solve the element.As shown in Figure 15, the bistatic RCS of GSST-EPA agrees very well with that of EFIE-PMCHWT.The comparison of computational results is given in Table 1.Because the number of unknowns on the equivalence surfaces is much less than that on the element surfaces, GSST-EPA can reduce the solution time and memory requirement significantly.It is well known that EFIE-PMCHWT is accurate but suffers from poor spectrum property.By transferring the unknowns on the element surfaces to the unknowns on the equivalence surfaces, GSST-EPA with CRM can improve the spectrum property of EFIE-PMCHWT significantly.
Finally, to investigate the capability of GSST-EPA with CRM to solve large array structures, the scattering of the 15 × 15 bowtie antenna array is solved.The parameters of the array are the same as that in the above example, except for the number of elements.The total number of unknowns of the equivalence surfaces is 41,850.The MLFMA is applied to accelerate the calculation of the coupling between nonadjacent elements.The bistatic RCS of HH polarization is given in Figure 16.The total solution time is 1972.6 s, and memory consumption is 3.9 GB.GSST-EPA converges to 10 −3 using 74 iteration steps.It can be seen that GSST-EPA still exhibits fast convergence as the increment of the element number.

Conclusions
In this paper, GSST-EPA with the CRM method is introduced to solve the array structures with connected elements.Firstly, the array structures are decomposed into subdomains by CRM.Then, the unknowns on the subdomains are transferred to the unknowns on the equivalence surfaces by GSST-EPA.Moreover, the meshes of subdomains can be the same with each other, which makes the scattering matrix solved and stored only once.The shift invariant of translation matrices is also used according to the characteristic of array structures.Using this strategy, the number of unknowns can be reduced and the conditioning of the final matrix equation is significantly better than the traditional PMCHWT and EFIE-PMCHWT.Therefore, the improved

Figure 3 :
Figure 3: The dielectric object decomposed into subdomains and each subdomain enclosed by an equivalence surface.

Figure 5 :
Figure 5: The translation operator of two nonadjacent subdomains.

Figure 6 :
Figure 6: The translation operator of two adjacent subdomains.

Figure 8 :Figure 9 :
Figure 8: (a) The model of a dielectric sphere.(b) Subdomain partitions of the dielectric sphere.

Figure 13 :
Figure 13: The convergence of GSST-EPA for the dielectric slab.

Table 1 :
The computational results of two methods.