AIMx: An Extended Adaptive Integral Method for the Fast Electromagnetic Modeling of Complex Structures
IIEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 1
AIMx: An Extended Adaptive Integral Methodfor the Fast Electromagnetic Modelingof Complex Structures
Shashwat Sharma,
Student Member, IEEE, and Piero Triverio,
Senior Member, IEEE
Abstract —Surface integral equation (SIE) methods are of greatinterest for the efficient electromagnetic modeling of variousdevices, from integrated circuits to antenna arrays. Existingacceleration algorithms for SIEs, such as the adaptive integralmethod (AIM), enable the fast approximation of interactions be-tween well-separated mesh elements. Nearby interactions involvethe singularity of the kernel, and must instead be computedaccurately with direct integration at each frequency of interest,which can be computationally expensive. In this work, a novelalgorithm is proposed for reducing the cost-per-frequency as-sociated with near-region computations for both homogeneousand layered background media. In the proposed extended AIM(AIMx), the SIE operators are decomposed into a frequency-independent term, which contains the singularity of the kernel,and a frequency-dependent term, which is a smooth function.The expensive near-region computations are only required forthe frequency-independent term, and can be reused at eachfrequency point, leading to significantly faster frequency sweeps.The frequency-dependent term is accurately captured via theAIM even in the near region, as confirmed through erroranalysis. The accuracy and efficiency of the proposed method aredemonstrated through numerical examples drawn from severalapplications, and CPU times are significantly reduced by factorsranging from three to 16.
Index Terms —Electromagnetic modeling, surface integralequations, adaptive integral method, fast frequency sweep.
I. I
NTRODUCTION F ULL-wave electromagnetic modeling is a crucial step inmany design workflows, ranging from the characterizationof on-chip interconnect networks, to the prediction and synthe-sis of antenna radiation profiles. The growing complexity ofmodern electronic and communication devices places a com-mensurate demand on full-wave simulation tools, especiallysince device characterization is often required over a widefrequency band.Volumetric methods such as the finite element method(FEM) [1] and volume integral equations (VIEs) [2]–[4] haveproved to be robust, but require a volumetric discretizationof conductive and dielectric objects. The FEM additionallyrequires a discretization of the space in between objects. Theserequirements lead to systems of equations with a large number
Manuscript received . . . ; revised . . . .S. Sharma and P. Triverio are with the Edward S. Rogers Sr. De-partment of Electrical & Computer Engineering, University of Toronto,Toronto, ON, M5S 3G4 Canada, e-mails: [email protected],[email protected] work was supported by the Natural Sciences and Engineering Re-search Council of Canada (Collaborative Research and Development Grantsprogram), by Advanced Micro Devices, and by CMC Microsystems. of unknowns. In contrast, surface integral equation (SIE)-based techniques require only a discretization of the surfaceof conductive and dielectric objects [5]–[10]. However, SIEmethods lead to a dense system matrix whose solution is com-putationlly prohibitive for large problems. Techniques basedon SIEs heavily rely on various acceleration algorithms [11]–[15], where the system of equations is solved iteratively. Theproducts between a system matrix A and a vector x arecomputed at each iteration as Ax = A NR x + A FR x . (1)where A NR and A FR contain entries of A corresponding tonearby (“near-region”) and well-separated (“far-region”) meshelements, respectively. Both A NR and A FR depend on theGreen’s function of the surrounding medium. The entries of A NR are dominant, and must be computed accurately withdirect integration. Far-region interactions in A FR are weaker,therefore the product A FR x can be computed in a fast butapproximate way, without having to assemble A FR .The fast multipole method (FMM) [11], [12] and its multi-level and mixed-form extensions [13], [14] are popular acceler-ation algorithms. These techniques were developed specificallyfor the homogeneous medium Green’s function. However,many practical applications involve objects embedded in strati-fied dielectric media, which may require the more complicatedmultilayer Green’s function (MGF) [16]. The extension ofFMM-based techniques to handle the MGF is not straight-forward, and requires either multiple complex exponentialterms [17] or numerical integration in the complex plane [18].Another class of acceleration algorithms arises from the useof an auxiliary coarse grid to compute A FR x quickly. Theadaptive integral method (AIM) [15] involves superimposinga regular grid on the structure, and projecting currents andcharges from the original mesh onto the auxiliary grid. Thetranslation invariance of the Green’s function is then exploitedby using Fast Fourier Transforms (FFTs) to accelerate thecomputation of A FR x . An advantage of the AIM is its kernel-independence, which makes it relatively easy to incorporatethe MGF [19]–[22]. The AIM has inspired several relatedmethods, such as the precorrected FFT (pFFT) [23] method,the sparse-matrix/canonical grid (SMCG) [24] method, andseveral others [25]–[28].In existing acceleration methods, generating A NR withdirect integration may still have a significant cost, particularlyfor densely-packed structures with fine features, which oftenrequire an extremely fine local mesh. This cost is further a r X i v : . [ c s . C E ] S e p EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 2 increased in the presence of layered media, since the MGFis significantly more expensive to compute than the homoge-neous Green’s function. Moreover, since the Green’s functionis singular when the source and observation points coincide,the entries of A NR must be computed carefully with advancedsingularity extraction or cancellation routines [29], which fur-ther add to the computational cost. Finally, the computation of A NR must be repeated at each frequency of interest. Methodsto interpolate simulation results across frequency have beenproposed [30]–[34], but still require discrete simulation dataat multiple frequencies to gather sufficient sample points. Thisis particularly onerous for structures whose terminal responsevaries sharply with frequency, for example due to resonances.In this contribution, we propose a novel acceleration algo-rithm based on the AIM, which requires computing A NR onlyonce for the entire frequency sweep. In the proposed extendedadaptive integral method (AIMx), the kernel of each integraloperator is carefully decomposed into a frequency-independentterm which contains the singularity of the kernel, and afrequency-dependent term which is smooth and nonsingulareverywhere in space. This leads to a convenient decompositionof the associated system matrix into frequency-independentand frequency-dependent parts. Through analytical and nu-merical error analyses, we show that the frequency-dependentpart of the system matrix can be accurately interpolatedwith polynomials even in the near region. Therefore, matrix-vector products involving the frequency-dependent part of eachmatrix operator are accelerated with the use of FFTs. Thefrequency- independent part of the matrix is expressed as in(1), where the near-region entries are computed accuratelywith direct integration. As a consequence, the key advantageof the proposed method is that direct integration needs to beperformed only once at the beginning of a frequency sweep,and then can be reused at each frequency point. We describein detail the proposed matrix decomposition for both thesingle- and double-layer potential matrix operators [35], forboth homogeneous and layered background media. This leadsto significantly faster frequency sweeps for a wide range ofproblems requiring full-wave electromagnetic characterization.The paper is organized as follows: Section II summarizesthe SIE formulation considered and the conventional AIM.Section III provides a detailed description of the proposedmethod, and Section IV shows analytical and numerical erroranalyses to support the hypothesis that frequency-dependentterms can be interpolated accurately in the near region. InSection V, we demonstrate the accuracy and efficiency of theproposed method compared to the conventional AIM for aseries of numerical examples, including canonical geometries,on-chip devices and large antenna arrays. A discussion of theadvantages of the proposed method is given in Section VI,followed by concluding remarks in Section VII.II. B ACKGROUND
Throughout this work, we consider time-harmonic fields. Atime dependence of e jωt is assumed and suppressed. Vectorquantities are written with an overhead arrow, for example (cid:126)a ( (cid:126)r ) . Primed coordinates represent source points, while un-primed coordinates represent observation points. Matrices and column vectors are written in bold letters, such as A and x ,while dyadic quantities are denoted with double overbars, asin A . A. Formulation
We consider the problem of electromagnetic scattering froma perfect electric conductor (PEC) object to illustrate theproposed method, though it is also applicable to penetrableobjects. The surface of the PEC object is denoted as S ,with outward unit normal ˆ n . The conductor is embedded infree space, denoted by V . The layered medium case will beconsidered in Section III-C. The unknown induced electricsurface current density (cid:126)J ( (cid:126)r ) on S due to an incident field, [ (cid:126)E inc ( (cid:126)r ) , (cid:126)H inc ( (cid:126)r )] , (cid:126)r ∈ V , satisfies the electric field integralequation (EFIE) [29], − jωµ ˆ n × L (cid:104) (cid:126)J ( (cid:126)r (cid:48) ) (cid:105) ( (cid:126)r ) = ˆ n × (cid:126)E inc ( (cid:126)r ) , (2)where the single-layer potential operator [35] L is definedas [29] L (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r )= ˆ S (cid:18) (cid:126)X ( (cid:126)r (cid:48) ) + 1 k ∇ (cid:48) ∇ (cid:48) · (cid:126)X ( (cid:126)r (cid:48) ) (cid:19) G ( k , (cid:126)r, (cid:126)r (cid:48) ) d S . (3)In (2) and (3), ω is the cyclical frequency, µ and k are thepermeability and wave number of free space, respectively, andthe kernel G ( k , (cid:126)r, (cid:126)r (cid:48) ) is the free space Green’s function, G ( k , (cid:126)r, (cid:126)r (cid:48) ) = e − jk r πr , (4)where r = | (cid:126)r − (cid:126)r (cid:48) | .Next, a triangular mesh is generated for S , and (2) isdiscretized with the method of moments (MoM). The electricsurface current density (cid:126)J ( (cid:126)r ) is expanded with Rao-Wilton-Glisson (RWG) basis functions [36] defined on pairs ofadjacent triangles, and (2) is tested with ˆ n × RWG basisfunctions to get the matrix equation − jωµ (cid:18) L ( A ) + 1 k L ( φ ) (cid:19) J = E inc , (5)where L ( A ) and L ( φ ) are, respectively, the vector and scalarpotential parts of the discretized ˆ n × L operator. Columnvectors J and E inc contain coefficients of the basis functionsassociated with (cid:126)J ( (cid:126)r ) and ˆ n × (cid:126)E inc ( (cid:126)r ) , respectively. Thehypersingular part of (3) is handled by transferring the gradientoperator ∇ (cid:48) to the testing function [29]. B. The Conventional AIM
In the AIM, a regular grid is superimposed on the boundingbox associated with S , and L is written as L = L NR + L FR , (6)where L represents either L ( A ) or L ( φ ) . Matrix L NR is sparseand contains the near-region entries of L , and is computed EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 3 accurately using direct integration. Matrix L FR contains far-region interactions, which are less dominant and can beapproximated as L FR ≈ W ( L ) H ( k ) P ( L ) − L P , (7)where matrix P ( L ) projects sources from the mesh to the AIMgrid, H ( k ) is the convolution matrix encoding the Green’sfunction, and W ( L ) interpolates the computed potentials backonto the mesh. Due to the Toeplitz nature of H ( k ) , matrix-vector products involving L FR are accelerated with FFTs [23].The AIM approximation is inaccurate in the near-region dueto the singularity of G ( k , (cid:126)r, (cid:126)r (cid:48) ) when (cid:126)r = (cid:126)r (cid:48) . Therefore, theprecorrection matrix L P is used to cancel out the contributionof W ( L ) H ( k ) P ( L ) in the near region [23], L P = W ( L ) H NR ( k ) P ( L ) , (8)where H NR ( k ) is the convolution matrix associated with thenear region.Although L NR is sparse, its assembly with direct integrationis still expensive for large problems, particularly when thestructure is densely packed and requires a fine mesh to resolveintricate features. This is especially cumbersome for layeredbackground media which require the expensive multilayerGreen’s function (MGF) [16], which is significantly moreexpensive to compute than its free space counterpart. Both L NR and L P must be generated at each frequency point, whichcan be a major bottleneck for wideband frequency sweeps.III. P ROPOSED M ETHOD
The conventional AIM has the following limitation: due tothe singularity of the Green’s function at r = 0 , the entriesof L NR cannot be approximated with interpolation, and mustinstead be computed with direct integration. The singularbehavior of G ( k , (cid:126)r, (cid:126)r (cid:48) ) may be investigated further througha Taylor expansion about r = 0 , G ( k , (cid:126)r, (cid:126)r (cid:48) ) = 14 πr + ( − jk )4 π + ( − jk ) π · r + ( − jk ) π · r + · · · . (9)The first term on the right-hand side of (9) is frequency-independent, and is the static Green’s function [35]. Theremaining terms on the right-hand side of (9) are frequency-dependent. Therefore, we may write G ( k , (cid:126)r, (cid:126)r (cid:48) ) = G s ( (cid:126)r, (cid:126)r (cid:48) ) + G d ( k , (cid:126)r, (cid:126)r (cid:48) ) , (10)where G s ( (cid:126)r, (cid:126)r (cid:48) ) = 14 πr (11)is frequency-independent, and G d ( k , (cid:126)r, (cid:126)r (cid:48) ) = G ( k , (cid:126)r, (cid:126)r (cid:48) ) − G s ( (cid:126)r, (cid:126)r (cid:48) ) = e − jk r − πr . (12)is frequency-dependent. Equations (9) and (10) reveal the keyobservations underlying the proposed method: • the singularity of the Green’s function is contained en-tirely in its frequency-independent part, G s ( (cid:126)r, (cid:126)r (cid:48) ) ; • since each frequency- dependent term on the right-handside of (9) is either constant, or decays to as r → , G d ( k , (cid:126)r, (cid:126)r (cid:48) ) is a smooth function for any finite r .This leads to our main hypothesis: G d ( k , (cid:126)r, (cid:126)r (cid:48) ) is amenableto approximation by interpolation, even in the near regionwhere r is small. Therefore, L NR may also be decomposedinto frequency-independent and frequency-dependent parts.Then, only the frequency-independent part needs to be com-puted accurately via direct integration, which does not have tobe repeated at each frequency point. Instead, the frequency-dependent part is computed in an accelerated manner with theAIM even in the near region , leading to the proposed extendedadaptive integral method (AIMx). A. AIMx for the Single-Layer Potential Operator
Using (10) in (3) allows decomposing the single-layerpotential operator L as L (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) = L s (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( (cid:126)r ) + L d (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) , (13)where the frequency-independent part is L s (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( (cid:126)r )= ˆ S (cid:18) (cid:126)X ( (cid:126)r (cid:48) ) + 1 k ∇ (cid:48) ∇ (cid:48) · (cid:126)X ( (cid:126)r (cid:48) ) (cid:19) G s ( (cid:126)r, (cid:126)r (cid:48) ) d S , (14)and the frequency-dependent part is L d (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r )= ˆ S (cid:18) (cid:126)X ( (cid:126)r (cid:48) ) + 1 k ∇ (cid:48) ∇ (cid:48) · (cid:126)X ( (cid:126)r (cid:48) ) (cid:19) G d ( k , (cid:126)r, (cid:126)r (cid:48) ) d S . (15)Correspondingly, L may be decomposed as L ( k ) = L s + L d ( k ) , (16)where L s and L d ( k ) represent either the vector or scalarpotential parts of the discretized L s and L d operators, respec-tively. The goal is to avoid the explicit assembly of L d , andto instead directly compute the product of L d with a vectorwith the AIM.To this end, only the frequency-independent matrix L s issplit into near- and far-region contributions as per (6), L ( k ) = L s , NR + L s , FR + L d ( k ) . (17)In (17), the singularity of G ( k , (cid:126)r, (cid:126)r (cid:48) ) is contained entirelyin the sparse, frequency-independent matrix L s , NR , whichis computed accurately with direct integration. The AIM isemployed to express the remaining matrices, L s , FR + L d ( k ) ≈ − L s , P + W ( L ) H s P ( L ) + W ( L ) H d ( k ) P ( L ) , (18)where L s , P = W ( L ) H s , NR P ( L ) , (19) EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 4 is the near-region precorrection matrix associated to L s , NR ,and H s and H d ( k ) are convolution matrices depending on G s ( (cid:126)r, (cid:126)r (cid:48) ) and G d ( k , (cid:126)r, (cid:126)r (cid:48) ) , respectively. Matrix H s , NR isthe convolution matrix for near-region precorrections, andencodes G s ( (cid:126)r, (cid:126)r (cid:48) ) . It should be noted that precorrections arerequired only for L s , NR , because both near- and far-regioncontributions of L d ( k ) are approximated with the AIM. Since H s + H d ( k ) = H ( k ) , (18) can be simplified further bycombining its last two terms, L ( k ) ≈ L s , NR − L s , P + W ( L ) H ( k ) P ( L ) . (20)Equation (20) is the proposed AIMx matrix decomposition forthe single-layer potential operator. Finally, using (20) in (5)yields the AIMx-accelerated EFIE for PECs, − jωµ (cid:104) L ( A )s , NR − L ( A )s , P + W ( A ) H ( k ) P ( A ) + 1 k (cid:16) L ( φ )s , NR − L ( φ )s , P + W ( φ ) H ( k ) P ( φ ) (cid:17)(cid:21) J = E inc , (21)where the superscript ( A ) or ( φ ) indicates that the associatedmatrix corresponds to the vector or scalar potential part of L ,respectively.In (21), all frequency-dependent effects are approximatedwith interpolation onto the auxiliary grid, and the associatedmatrix-vector products are accelerated with FFTs. The keyadvantage of the proposed method is that the entries of L ( A )s , NR , L ( φ )s , NR , L ( A )s , P and L ( φ )s , P are frequency-independent, andneed only be computed once, enabling significantly fasterfrequency sweeps compared to the conventional AIM. Animportant property of the proposed method is its simplicityof implementation. Compared to an existing implementationof the conventional AIM, the proposed method requires onlytwo changes:1) the near-region and precorrection matrices, L s , NR and L s , P , must be assembled using the kernel G s ( (cid:126)r, (cid:126)r (cid:48) ) rather than G ( k , (cid:126)r, (cid:126)r (cid:48) ) ;2) a minor modification to H must be made, which isdescribed in Section III-A1The fact that direct integration is only required for thefrequency-dependent part of (10) was also exploited in [37],where the Taylor expansion (9) is explicitly used to obtain aseries of matrices to accurately handle very low frequencies.However, unlike [37], the proposed method does not explicitlyuse a series expansion, and therefore avoids the need for anincreasingly large number of matrix terms at high frequencies.Furthermore, the proposed method is applicable to both thesingle- and double-layer potential operators [35] for both thefree space and the multilayer Green’s function. The Taylorexpansion of the Green’s function has also been used to speedup near-region computations in [38] for free space, and [39],[40] for layered media. However, unlike these methods, theproposed technique does not explicitly utilize a Taylor expan-sion, and allows the direct integration step to be performedindependently of the frequency.
1) Modification to the Convolution Matrix H : Assumingthat the AIM grid contains N g grid points, the ( m, n ) entryof H is H ( m,n ) = G ( k , (cid:126)r m , (cid:126)r n ) , (22)where m, n = 1 , . . . N g , and (cid:126)r i is the i th grid point. Due tothe regularity of the AIM grid and the translation-invarianceof G ( k , (cid:126)r, (cid:126)r (cid:48) ) , H is a three-level Toeplitz matrix [23].Therefore, only its first row contains unique entries which mustbe computed and stored. When m = n , the source and obser-vation grid points coincide, i.e. r = 0 . The term H ( m,m ) is notwell defined, since G ( k , (cid:126)r, (cid:126)r (cid:48) ) → ∞ as r → . However, the r = 0 case necessarily resides within the near region, and isaccounted for with direct integration using appropriate singu-larity subtraction or cancellation techniques [29]. Therefore,in the conventional AIM, one can set H ( m,m ) to an arbitraryvalue, so long as H ( m,m )NR is set to the same value. This ensuresthat H ( m,m ) is cancelled out when near-region precorrectionsare applied.In the proposed AIMx, the H ( m,m ) term must be handledmore carefully. The contribution of G s ( (cid:126)r, (cid:126)r (cid:48) ) in the nearregion, particularly when r = 0 , is accounted for by L s , NR via direct integration, with the extraction or cancellation ofthe singularity [29]. Therefore, H ( m,m )s can be set to anarbitrary value, so long as H ( m,m )s , NR is set to the same value.For simplicity, we set H ( m,m )s = H ( m,m )s , NR = 0 . (23)The term H ( m,m )d contains the contribution of G d ( k , (cid:126)r, (cid:126)r (cid:48) ) for r = 0 , which must still be accounted for and can beobtained by inspection of (9). The first term on the right-handside of (9) is G s ( (cid:126)r, (cid:126)r (cid:48) ) , while the remaining terms represent G d ( k , (cid:126)r, (cid:126)r (cid:48) ) , G d ( k , (cid:126)r, (cid:126)r (cid:48) ) = ( − jk )4 π + ( − jk ) π · r + ( − jk ) π · r + · · · . (24)When r = 0 , all but the first term on the right-hand side of (24)are . Therefore, we let H ( m,m )d = G d ( k , (cid:126)r m , (cid:126)r m ) = ( − jk )4 π . (25)Finally, since H ( m,m ) = H ( m,m )s + H ( m,m )d , using (23)and (25) we have H ( m,m ) = H ( m,m )d = ( − jk )4 π , (26)which is the simple modification that must be made to H compared to the conventional AIM. This modification ensuresthat there is no undefined behavior when r = 0 , and that H correctly accounts for frequency-dependent effects for allvalues of r , including r = 0 .
2) Preconditioning:
The fact that only frequency-independent near-region interactions are precomputed withthe MoM has important implications for constructing apreconditioner to solve (21) iteratively. Common choicesof preconditioner include diagonal, near-region and block-diagonal entries of L [41], [42]. The key challenge in EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 5 preconditioner assembly for the proposed method is that L d ( k ) is not assembled explicitly, and only the entries of L s are available prior to the iterative solution of (21).In the case of diagonal preconditioning, only the self-termsof L are involved, which are dominated by the contribution of G s ( (cid:126)r, (cid:126)r (cid:48) ) . Therefore, it is reasonable to expect that in the pro-posed method, using the diagonal of L s as the preconditionershould be just as effective as using the diagonal of L , as in theconventional AIM. This is indeed the case, as demonstratedthrough several numerical examples in Section V.In the case of near-region or block-diagonal precondition-ing, using the entries of L s in the proposed AIMx may notbe as effective as using those of L , as in the conventionalAIM. This is because the off-diagonal entries of L becomeincreasingly influenced by G d ( k , (cid:126)r, (cid:126)r (cid:48) ) as r increases. Ifoff-diagonal entries of L are required for preconditioning,one may use FFTs via (20) to explicitly generate selectedentries of L d ( k ) to augment the corresponding entries of L s .This procedure may be performed at each frequency point, orat every few frequency points. The associated computationalcost is similar to that of the near-region precorrection phaseof the conventional AIM, which may be significant but isnot typically a bottleneck. Crucially, the direct integrationprocedure is still avoided at each new frequency point. B. AIMx for the Double-Layer Potential Operator
When S is a closed surface, one may also express (cid:126)J ( (cid:126)r ) viathe magnetic field integral equation (MFIE) [29], − ˆ n × ˆ n × K (cid:104) (cid:126)J ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) = ˆ n × ˆ n × (cid:126)H inc ( (cid:126)r ) , (27)where, in free space, the double-layer potential operator [35] K is defined as [43] K (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) = ˆ S ∇ G ( k , (cid:126)r, (cid:126)r (cid:48) ) × (cid:126)X ( (cid:126)r (cid:48) ) d S , (28)where ∇ G ( k , (cid:126)r, (cid:126)r (cid:48) ) = − ˆ r (1 + jk r ) e − jk r πr , (29)and ˆ r = ( (cid:126)r − (cid:126)r (cid:48) ) /r . For closed PEC objects, it is oftendesirable to formulate an integral equation by adding (2) toa scaled version of (27), to obtain the combined field integralequation (CFIE) [29], − jωµ ˆ n ×L (cid:104) (cid:126)J ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) − η ˆ n × ˆ n ×K (cid:104) (cid:126)J ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r )= ˆ n × (cid:126)E inc ( (cid:126)r ) + η ˆ n × ˆ n × (cid:126)H inc ( (cid:126)r ) , (30)where η is the wave impedance of free space. Discretizingand testing (30) as before yields its discrete counterpart, [ − jωµ L − η K ] J = E inc + η H inc , (31)where K is the discretized ˆ n × ˆ n × K operator, and columnvector H inc contains coefficients of the basis functions asso-ciated with ˆ n × (cid:126)H inc ( (cid:126)r ) . To develop the proposed method for the K operator, weexpand ∇ G ( k , (cid:126)r, (cid:126)r (cid:48) ) in a Taylor series about r = 0 , ∇ G ( k , (cid:126)r, (cid:126)r (cid:48) )= − ˆ r π (cid:32) r − ( − jk ) − − jk ) r + · · · (cid:33) . (32)As in the case of G ( k , (cid:126)r, (cid:126)r (cid:48) ) , the singularity of ∇ G ( k , (cid:126)r, (cid:126)r (cid:48) ) is contained entirely in the first term on the right-handside of (32), which is frequency-independent. The remainingfrequency-dependent terms are either constant or decay to as r → , and are smooth functions for small r . Following thesame approach devised in Section III-A for L , we may express K as K (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) = K s (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( (cid:126)r ) + K d (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) , (33)where the frequency-independent part is K s (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( (cid:126)r ) = ˆ S ∇ G s ( (cid:126)r, (cid:126)r (cid:48) ) × (cid:126)X ( (cid:126)r (cid:48) ) d S , (34)and the frequency-dependent part is K d (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( k , (cid:126)r ) = ˆ S ∇ G d ( k , (cid:126)r, (cid:126)r (cid:48) ) × (cid:126)X ( (cid:126)r (cid:48) ) d S , (35)In (34) and (35), ∇ G s ( (cid:126)r, (cid:126)r (cid:48) ) = − ˆ r πr , (36)and ∇ G d ( k , (cid:126)r, (cid:126)r (cid:48) ) = ∇ G ( k , (cid:126)r, (cid:126)r (cid:48) ) − ∇ G s ( (cid:126)r, (cid:126)r (cid:48) ) . (37)Therefore, K may be decomposed as K ( k ) = K s + K d ( k ) , (38)where K s and K d ( k ) are the discretized K s and K d operators,respectively.Following the approach in Section III-A, K may be approx-imated as K ≈ K s , NR − K s , P + W ( K ) H ∇ P ( K ) , (39)where K s , NR and K s , P are the near-region and precorrectionmatrices, respectively, which are frequency-independent andinvolve the kernel ∇ G s ( (cid:126)r, (cid:126)r (cid:48) ) , and H ∇ is the convolutionmatrix encoding ∇ G ( k , (cid:126)r, (cid:126)r (cid:48) ) . Matrices W ( K ) and P ( K ) arethe interpolation and projection matrices associated with K ,respectively. In keeping with the discussion in Section III-A1,the self-term of H ∇ should be set as H (0 , ∇ = ∇ G d ( k , (cid:126)r , (cid:126)r ) = ( − jk ) π · , (40)which is the second term of the expansion (32). The singularityof ∇ G s ( (cid:126)r, (cid:126)r (cid:48) ) at r = 0 is accounted for with direct integrationvia K s , NR . EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 6
Finally, using (20) and (39) in (31), we obtain the AIMx-accelerated CFIE for PECs, − jωµ (cid:20) (cid:16) L ( A )s , NR − L ( A )s , P + W ( A ) HP ( A ) (cid:17) + 1 k (cid:16) L ( φ )s , NR − L ( φ )s , P + W ( φ ) HP ( φ ) (cid:17) − jk (cid:16) K s , NR − K s , P + W ( K ) H ∇ P ( K ) (cid:17)(cid:21) J = E inc + η H inc . (41)The comments regarding preconditioners in Section III-A2apply to (41) as well. C. AIMx for Layered Background Media
For layered background media, we consider the multilayerGreen’s function (MGF) as described in formulation C of [44],though the proposed method may be extended to other formu-lations as well. The definition of L becomes [43] L (cid:104) (cid:126)X ( (cid:126)r (cid:48) ) (cid:105) ( (cid:126)r ) = ˆ S G A ( k , (cid:126)r, (cid:126)r (cid:48) ) · (cid:126)X ( (cid:126)r (cid:48) ) d S + 1 k ˆ S ∇ (cid:48) ∇ (cid:48) · (cid:126)X ( (cid:126)r (cid:48) ) G φ ( k , (cid:126)r, (cid:126)r (cid:48) ) d S , (42)where G A ( k , (cid:126)r, (cid:126)r (cid:48) ) and G φ ( k , (cid:126)r, (cid:126)r (cid:48) ) are the dyadic andscalar potential Green’s functions of the background medium,respectively. The term G φ ( k , (cid:126)r, (cid:126)r (cid:48) ) and the non-zero com-ponents of G A ( k , (cid:126)r, (cid:126)r (cid:48) ) are expressed in terms of semi-infinite integrals in the complex plane, which are knownas “Sommerfeld integrals” [16], [45]. These integrals areexpensive to compute, even when approximations such as thediscrete complex image method (DCIM) [46] are employed.As a result, the near-region computations involving directintegration are particularly time consuming for layered media.In the following, it is assumed that G ( k , (cid:126)r, (cid:126)r (cid:48) ) rep-resents either G φ ( k , (cid:126)r, (cid:126)r (cid:48) ) , or a non-zero component of G A ( k , (cid:126)r, (cid:126)r (cid:48) ) . Regardless of the method used to compute G ( k , (cid:126)r, (cid:126)r (cid:48) ) , it is common practice to extract the quasistaticcontributions in the spectral domain, which are then addedback analytically in spatial domain [45]–[47]. This ensuresthat the integrands involved in the computation of G ( k , (cid:126)r, (cid:126)r (cid:48) ) decay sufficiently fast [47]. Extracting quasistatic terms alsoprovides a simple way to compute the frequency-independentpart of G ( k , (cid:126)r, (cid:126)r (cid:48) ) , as discussed below. Therefore, we maywrite G ( k , (cid:126)r, (cid:126)r (cid:48) ) = G q ( k , (cid:126)r, (cid:126)r (cid:48) ) + G r ( k , (cid:126)r, (cid:126)r (cid:48) ) , (43)where G q ( k , (cid:126)r, (cid:126)r (cid:48) ) is the quasistatic part of G ( k , (cid:126)r, (cid:126)r (cid:48) ) , and G r ( k , (cid:126)r, (cid:126)r (cid:48) ) is the remainder. In general, G q ( k , (cid:126)r, (cid:126)r (cid:48) ) maybe expressed as [47] G q ( k , (cid:126)r, (cid:126)r (cid:48) ) = N s (cid:88) n =1 G s n ( (cid:126)r, (cid:126)r (cid:48) ) + N q (cid:88) n =1 G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) , (44)where N s , N q and the forms of G s n ( (cid:126)r, (cid:126)r (cid:48) ) and G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) all depend on the location of source and observationpoints [47]. In (44), G q ( k , (cid:126)r, (cid:126)r (cid:48) ) contains some terms, G s n ( (cid:126)r, (cid:126)r (cid:48) ) , that are already frequency-independent. The otherterms, G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) , depend on frequency, and have theform [47], G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) = C n π e − jk i α n α n , (45)where C n is a constant which depends on the material param-eters of the layered medium [47], k i is the wave number inthe layer in which the source resides, ρ = (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) , (46)and α n = (cid:112) ρ + γ n . (47)Term γ n is a linear function of z , z (cid:48) and the elevation of sourceand observation layers [47]. Equation (45) has a form similarto (4), and the Taylor expansion of G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) about α = 0 reads G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) = C n (cid:20) πα n + ( − jk i )4 π + ( − jk i ) π · α n + ( − jk i ) π · α n + · · · (cid:35) , (48)where again the singularity is contained in the first term, whichis frequency-independent. We may decompose G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) into frequency-independent and frequency-dependent parts as G q n ( k , (cid:126)r, (cid:126)r (cid:48) ) = C n πα n (cid:124) (cid:123)(cid:122) (cid:125) G qs n ( (cid:126)r,(cid:126)r (cid:48) )+ C n π (cid:18) e − jk i α n α n − α n (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) G qd n ( k ,(cid:126)r,(cid:126)r (cid:48) ) . (49)Using (49) and (44) in (43), we get G ( k , (cid:126)r, (cid:126)r (cid:48) ) = N s (cid:88) n =1 G s n ( (cid:126)r, (cid:126)r (cid:48) ) + N q (cid:88) n =1 G qs n ( (cid:126)r, (cid:126)r (cid:48) )+ N q (cid:88) n =1 G qd n ( k , (cid:126)r, (cid:126)r (cid:48) ) + G r ( k , (cid:126)r, (cid:126)r (cid:48) ) , (50)where the first two terms on the right-hand side are frequency-independent and can be denoted as G s ( (cid:126)r, (cid:126)r (cid:48) ) , while thelast two terms are frequency-dependent and are denoted as G d ( k , (cid:126)r, (cid:126)r (cid:48) ) . Therefore, (50) may be simplified as G ( k , (cid:126)r, (cid:126)r (cid:48) ) = G s ( (cid:126)r, (cid:126)r (cid:48) ) + G d ( k , (cid:126)r, (cid:126)r (cid:48) ) , (51)which resembles the free space case in (10). This expressionholds for both G φ ( k , (cid:126)r, (cid:126)r (cid:48) ) and any non-zero component of G A ( k , (cid:126)r, (cid:126)r (cid:48) ) . Equation (51) indicates that the main ideas ofthe proposed method also apply to layered media: • the singular behavior of G ( k , (cid:126)r, (cid:126)r (cid:48) ) is entirely containedin the frequency-independent term G s ( (cid:126)r, (cid:126)r (cid:48) ) ; • the frequency-dependent term G d ( k , (cid:126)r, (cid:126)r (cid:48) ) is smooth,and decays to as r → [45], [47], therefore can beapproximated by interpolation even in the near region.When the quasistatic part of the MGF is expressed as performulation C in [44], C n may depend on the complex permit-tivity of the associated layer, which is frequency-dependent if EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 7 the layer is lossy. As a result, for lossy layers, G qs n ( (cid:126)r, (cid:126)r (cid:48) ) mayno longer be frequency-independent. A complete treatment ofthe lossy case is beyond the scope of this paper, and will be thesubject of a future publication. In this work, we assume that thelayers are either lossless, or that their losses are small enoughthat the imaginary part of the complex permittivity can beneglected when (45) is evaluated in the near region. As shownin the numerical example in Section V-E, this approximationstill yields accurate results for realistic structures of practicalimportance.The AIMx concept may also be extended to the curlof G A ( k , (cid:126)r, (cid:126)r (cid:48) ) in order to handle the K operator forlayered media [43]. The proposed AIMx-based EFIE, (21),and CFIE, (41), can therefore be applied for PEC objectsembedded in layered media. An important advantage of theproposed method for layered media is that the near-regionmatrix entries are extremely simple to compute, by takingthe static limit of G q ( k , (cid:126)r, (cid:126)r (cid:48) ) . These entries do not requirecomputing Sommerfeld integrals, nor do they require applyingapproximations such as the DCIM. Therefore, the near-regioncomputations in the proposed method are significantly moreefficient, robust, and easy to implement, compared to theconventional AIM. IV. E RROR A NALYSIS
The proposed method rests on the hypothesis that thefrequency-dependent part of the kernels of L and K canbe approximated with interpolation in both the near and farregions. In this section, this hypothesis is validated for eachtype of kernel considered, through analytical and numericalerror analysis.We consider the moment-matching version of the AIM [48],which relies on the interpolation of the kernel with poly-nomials [23]. We use Lagrange polynomials due to theirstraightforward implementation and the availability of an errorbound. Without loss of generality, we can assume that thesource is at the origin, so that (cid:126)r (cid:48) = 0 , and r = | (cid:126)r | . In addition,we normalize any physical length, say r , as θ r (cid:44) k r . A. Free Space
We first consider the free space case, and assume that K ( θ r ) can represent either G d ( θ r ) or ∇ G d ( θ r ) . Consider aregular AIM grid with a stencil which has n + 1 points in eachdirection. Given n + 1 samples, ( r , K ( θ r )) , ( r , K ( θ r )) . . . ( r n , K ( θ r n )) , (52)the Lagrange polynomial approximation of K ( θ r ) can bewritten as [49] (cid:101) K n ( θ r ) = n (cid:88) q =0 K ( θ r q ) m (cid:54) = q (cid:89) ≤ m ≤ n r − r m r q − r m . (53)The order n interpolation error is ∆ K n ( θ r ) = (cid:12)(cid:12)(cid:12) K ( θ r ) − (cid:101) K n ( θ r ) (cid:12)(cid:12)(cid:12) . (54)The goal is to show that ∆ K n ( θ r ) is bounded even when θ r → , and is controllable via grid refinement. Following the approach in [50], ∆ K n ( θ r ) can be boundas [49] ∆ K n ( θ r ) k ≤ ( θ r n − θ r ) n +1 ( n + 1)! max θ r ≤ θ ξ ≤ θ rn (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) d ( n +1) dθ ( n +1) ξ [ K ( θ ξ )] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) B n ( θ r ) , (55)where θ r and θ r n are the stencil coordinates near-est to and farthest from the source point, respectively.When K ( θ r ) = G d ( θ r ) , d ( n +1) dθ ( n +1) ξ [ K ( θ ξ )] = ( − n +1 πθ n +2 ξ ( n + 1)! · (cid:34) − n +1 (cid:88) l =0 (cid:18) n + 1 l (cid:19) ( n + 1 − l )!( n + 1)! ( − jθ ξ ) l e − jθ ξ (cid:35) , (56)and when K ( θ r ) = ∇ G d ( θ r ) , d ( n +1) dθ ( n +1) ξ [ K ( θ ξ )] = − ˆ θ π n +1 (cid:88) l =0 (cid:18) n + 1 l (cid:19) ( − n − l ) ( n +1 − l ) θ n +3 − lξ · (cid:2) − δ l + ( − j ) l e − jθ ξ (1 + jθ ξ − l ) (cid:3) . (57)In (57), ˆ θ = k ˆ r , the notation ( a ) ( b ) indicates the risingfactorial (“Pochhammer function”) [49], and δ l = (cid:40) if l = 0 , otherwise. (58)When θ r = 0 , (56) and (57) cannot be evaluated numericallydue to their indeterminate forms. Instead, one can evaluate lim θ r → B n ( θ r ) by recursively applying the L’Hˆopital rule,or by using symbolic math software [51], to get lim θ r → B n ( θ r ) = ( θ r n − θ r ) n +1 ( n + 1)! (cid:18) n + 2 (cid:19) (59)for K ( θ r ) = G d ( θ r ) , and lim θ r → B n ( θ r ) = ( θ r n − θ r ) n +1 ( n + 1)! (cid:18) n + 3 (cid:19) (60)for K ( θ r ) = ∇ G d ( θ r ) . In both cases, lim θ r → B n ( θ r ) is aconstant, which confirms that the Lagrange interpolation errorfor both G d ( θ r ) and ∇ G d ( θ r ) is bounded even for θ r = 0 .Fig. 1 shows, over a wide range of electrical sizes, theresults of the error analysis for G d ( θ r ) , for n = 2 . Fig. 2shows the same for ∇ G d ( θ r ) . In both cases, the top panelshows the exact and interpolated curves, for the given samplepoints indicated with vertical black bars. It is assumed that samples are used per wavelength. The magenta dashed lineshows a typical choice for the size of the near region. Thebottom panel shows the computed numerical error, as well asthe analytical bound B n ( θ r ) . The green dashed line indicatesthe value of lim θ r → B n ( θ r ) obtained from (59) and (60).Fig. 1 and Fig. 2 clearly indicate that the error is boundedas θ r → , and is not significantly larger within the nearregion, as compared to the region just beyond. Also plottedin the bottom panel is the computed numerical error whenthe number of samples is increased to and points per EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 8 θ r / π | G d ( θ r ) | ExactInterpolatedInterpolation dataEdge of near-region θ r / π − − | ∆ K ( θ r ) k | B n ( θ r ), 10 pts/ λ lim θ r → B n ( θ r )Actual error, 10 pts/ λ Actual error, 20 pts/ λ Actual error, 30 pts/ λ Fig. 1: Top panel: Lagrange interpolation of G d ( θ r ) for n = 2 . Bottom panel: computed numerical error, and B n ( θ r ) computed via (56). θ r / π | ∇ G d ( θ r ) | ExactInterpolatedInterpolation dataEdge of near-region θ r / π − − | ∆ K ( θ r ) k | B n ( θ r ), 10 pts/ λ lim θ r → B n ( θ r )Actual error, 10 pts/ λ Actual error, 20 pts/ λ Actual error, 30 pts/ λ Fig. 2: Top panel: Lagrange interpolation of ∇ G d ( θ r ) for n = 2 . Bottom panel: computed numerical error, and B n ( θ r ) computed via (57).wavelength, confirming that the error is controllable via gridrefinement. Therefore, the frequency-dependent part of the freespace Green’s function and its gradient can be interpolatedwith polynomials for all θ r , including the near region. B. Layered Media
Deriving an analytical bound from (55) for layered mediais more challenging due to the many possible forms thatindividual dyadic components assume. Instead, we computethe bounds by evaluating the derivative in (55) with a finitedifference approximation for some representative examples.Table I shows the dielectric layer configuration used forthe error analysis in this section, and it is assumed that z = z (cid:48) = 0 . mm, where z = 0 is taken as the interface be-tween the lowermost layer and the PEC backing. Fig. 3 showsthe Lagrange interpolation and associated error where K ( θ r ) θ r / π | G d ( θ r ) | ExactInterpolatedInterpolation dataEdge of near-region θ r / π − − − | ∆ K ( θ r ) k | B n ( θ r ), 10 pts/ λ Actual error, 10 pts/ λ Actual error, 20 pts/ λ Actual error, 30 pts/ λ Fig. 3: Top panel: Lagrange interpolation of the frequency-dependent part of G A ( k , (cid:126)r, (cid:126)r (cid:48) ) , for the xx component, for n = 2 . Bottom panel: computed numerical error, and B n ( θ r ) computed via finite difference-based evaluation of (55).TABLE I: Dielectric layer configurations for the error analysisin Section IV-B and the numerical examples in Section V.Layers are non-magnetic, and h represents their height. Error Analysis Dipole Array Inductor Array(Section IV-B) (Section V-D) (Section V-E) ε r h (mm) ε r h (mm) ε r σ (S/m) h ( µ m)Air ∞ Air ∞ Air ∞ . . . . . . .
17 0 .
254 2 . . . .
04 5 .
77 3 . . . . . PEC ∞ PEC ∞ PEC ∞ is taken as the xx component of the frequency-dependentpart of G A ( k , (cid:126)r, (cid:126)r (cid:48) ) . As in the free space case, the errorremains bounded as θ r → , and is not significantly largerinside the near region than just outside. Furthermore, the errorcan be controlled via grid refinement. Therefore, the idea ofinterpolating dynamic contributions of the kernel in both nearand far regions extends to the modeling of layered media withthe MGF. V. R ESULTS
The accuracy and efficieny of the proposed method aredemonstrated in this section, through comparisons with theconventional AIM. All simulations were performed single-threaded on a 3 GHz Intel Xeon CPU. We used PETSc [52] forsparse matrix manipulation. The GMRES iterative solver [53]available through PETSc was used for solving (21) and (41).A relative residual norm of − was used as the GMRESconvergence tolerance. EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 9
TABLE II: Performance comparison for the numerical examples in Section V.
Sphere NASA Almond Reflectarray Dipole Array Inductor Array(Section V-A) (Section V-B) (Section V-C) (Section V-D) (Section V-E)AIM AIMx AIM AIMx AIM AIMx AIM AIMx AIM AIMxNumber of mesh elements ,
080 40 ,
194 254 ,
736 25 ,
042 42 , One-time static matrix-fill - - . sec - - . min - - . hr - - . min - - . hrPer-frequency matrix-fill* . sec - - . min - - . hr - - . min - - . hr - -Total CPU time . min . min . min . min . hr . hr . hr . hr . hr . hrOverall speed-up . × . × . × . × . × *averaged across all frequency points A. Sphere
We first consider a PEC sphere with diameter m in freespace. The sphere is meshed with , triangles and , edges, and an AIM grid with × × points along x , y and z , respectively, was used with n = 2 . A plane wave withthe electric field polarized along the y axis is incident on thesphere, and the EFIE system (21) was solved with a diagonalpreconditioner. The geometry and computed electric surfacecurrent density at MHz are shown in Fig. 4. The top panelof Fig. 5 shows the monostatic radar cross section (RCS)across a range of frequencies computed via the proposedmethod, compared to the conventional AIM and the analyticalMie series. Excellent agreement between the three methods isobserved across the entire frequency range. The bottom panelof Fig. 5 shows the bistatic RCS for the φ = 0 cut, for thelow ( MHz) and high (
MHz) ends of the consideredfrequency range. Again, the proposed method agrees perfectlywith the conventional AIM.The top panel of Fig. 6 shows that the proposed methodrequires the same number of iterations as the conventionalAIM except at two frequency points, which are near internalresonances. Therefore, the frequency-independent diagonalpreconditioner discussed in Section III-A2 is nearly as effec-tive as the full-wave diagonal preconditioner across the entirefrequency range. The bottom panel of Fig. 6 shows the CPUtime per frequency, and demonstrates the time saved by avoid-ing the near-region direct integration and precorrection steps ateach frequency. The time taken for near-region computationsin both AIMx and the conventional AIM is reported in Table II,and a total speed-up of . × was observed. B. NASA Almond
Next, we consider the NASA almond [54] in free space,meshed with , triangles and , edges. An AIM gridwith × × points was used with n = 2 . A plane wavewith the electric field polarized along the x axis is incidenton the PEC structure, and the CFIE system (41) was solvedwith a diagonal preconditioner. The geometry and computedelectric surface current density at GHz are shown in Fig. 7.The top pane of Fig. 8 shows the monostatic RCS across arange of frequencies computed via AIMx, compared to theconventional AIM. Excellent agreement is observed across theentire frequency range. The bottom panel of Fig. 8 showsthe bistatic RCS for the φ = 0 cut, at GHz and GHz.Again, the proposed method is in excellent agreement withthe conventional AIM. The top panel of Fig. 9 shows that the proposed methodrequires the same number of iterations to converge as theconventional AIM across the entire frequency range, whichreaffirms the effectiveness of the frequency-independent diag-onal preconditioner compared to the full-wave version. Thebottom panel of Fig. 9 shows the CPU time per frequency,and shows a significant reduction in the computational costcompared to the conventional AIM. The proposed methodprovided a total speed-up of . × , and CPU times are reportedin Table II. C. Planar Reflectarray
We next consider a more realistic test case, involving afour-layer reflectarray in free space with × planar PECJerusalem cross elements. Details regarding the geometry ofindividual elements can be found in [55], [56]. The array ismeshed with , triangles and , edges, and anAIM grid with × × points was used with n = 3 . Aplane wave with the electric field polarized along x impingeson the array. To robustly handle the intricate sub-wavelengthfeatures of the structure, the augmented EFIE (AEFIE) [7]was used, which involves the same matrix operators as theEFIE (5), but takes the charge density as an additional un-known. The constraint preconditioner [7] involving self-termsof L ( A )s , NR and L ( φ )s , NR was used. The geometry and computedelectric surface current density at GHz are shown in Fig. 10.The monostatic RCS for frequencies from to GHz isshown in the top panel of Fig. 11, which corresponds toelectrical sizes from . to . wavelengths. Additionally, thedirectivity for the φ = 0 cut is shown in the bottom panelof Fig. 11 for , and GHz. The RCS and directivitycomputed via the proposed method are both in excellentagreement with results from the conventional AIM.The top panel of Fig. 12 shows that once again, the precon-ditioner based on frequency-independent self-terms of L ( A )NR and L ( φ )NR is just as effective as the full-wave preconditionerproposed in [7] across the entire frequency range considered.Furthermore, the bottom panel of Fig. 12 shows the signifi-cantly improved efficieny of the proposed method comparedto the conventional AIM. At the lower end of the frequencyrange, fewer iterations are required in the iterative solve step,and the CPU time in the conventional AIM is dominated by thenear-region direct integration and precorrection steps. Instead,at higher frequencies, the number of iterations required toconverge increases, showing a commensurate increase in totalCPU time in the proposed method, as shown in the bottom EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 10
Fig. 4: Geometry and electric surface current density at
MHz for the sphere in Section V-A. f (MHz) − − − R C S ( d B s m ) MieConventional AIMProposed θ (degrees) − − B i s t a t i c R C S ( d B s m ) f = 10 MHz f = 10 MHz f = 750 MHz f = 750 MHz Conventional AIMProposed
Fig. 5: Accuracy validation for the sphere in Section V-A. Toppanel: monostatic RCS. Bottom panel: bistatic RCS at MHzand
MHz. f (MHz)200400600 I t e r a t i o n s Conventional AIMProposed f (MHz)102030 C P U t i m e ( s ) Conventional AIMProposed
Fig. 6: Performance comparison for the sphere in Section V-A.Top panel: iterations required for convergence. Bottom panel:CPU time per frequency.panel of Fig. 12. However, the overall cost of the conventionalAIM is still dominated by the near-region computations. TheCPU times for both methods are reported in Table II. It isimportant to note that besides the per-frequency time savings,the near-region direct integration step in the proposed method Fig. 7: Geometry and electric surface current density at GHzfor the NASA almond in Section V-B. f (GHz) − − − − R C S ( d B s m ) Conventional AIMProposed θ (degrees) − − B i s t a t i c R C S ( d B s m ) f = 1 GHz f = 1 GHz f = 10 GHz f = 10 GHz Conventional AIMProposed
Fig. 8: Accuracy validation for the NASA almond in Sec-tion V-B. Top panel: monostatic RCS. Bottom panel: bistaticRCS at GHz and GHz. f (GHz)406080100 I t e r a t i o n s Conventional AIMProposed f (GHz)0200400 C P U t i m e ( s ) Conventional AIMProposed
Fig. 9: Performance comparison for the NASA almond inSection V-B. Top panel: iterations required for convergence.Bottom panel: CPU time per frequency.is nearly × faster than in the conventional AIM. This isbecause the static Green’s function, / (4 πr ) , is used as thekernel for near-region computations in the proposed method,which is faster to compute than the exponential function in thefull-wave Green’s function (4), as required in the conventional EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 11
Fig. 10: Geometry and electric surface current density at GHz for the reflectarray in Section V-C.
16 18 20 22 24 f (GHz) − − R C S ( d B s m ) Conventional AIMProposed − − −
25 0 25 50 75 θ (degrees) − − D i r ec t i v i t y ( d B i ) f = 15 GHz f = 15 GHz f = 20 GHz f = 20 GHz f = 25 GHz f = 25 GHz Conventional AIMProposed
Fig. 11: Accuracy validation for the reflectarray in Sec-tion V-C. Top panel: monostatic RCS. Bottom panel: direc-tivity at GHz, GHz and GHz.
16 18 20 22 24 f (GHz)100200300 I t e r a t i o n s Conventional AIMProposed
16 18 20 22 24 f (GHz)10 C P U t i m e ( s ) Conventional AIMProposed
Fig. 12: Performance comparison for the reflectarray in Sec-tion V-C. Top panel: iterations required for convergence.Bottom panel: CPU time per frequency.AIM. The proposed method reduces the total simulation timefrom over days to just . hours, which corresponds to an overall speed-up of . × . D. Printed Dipole Array
The numerical examples considered so far involved struc-tures in free space; the following examples concern structuresembedded in layered media. First, we consider a double-sidedprinted dipole array consisting of planar PEC elements fedby a balanced transmission line network. The structure isshown in Fig. 13, and was inspired by [57], which includesdetails of the geometry. Each arm of the dipole elements isprinted on opposide sides of a dielectric substrate with relativepermittivity . , adjacent to a foam spacer backed by aground plane. The background medium configuration is shownin Table I. The structure is meshed with , trianglesand , edges, and an AIM grid with × × points was used with n = 3 . The AEFIE formulation withthe same constraint preconditioner as in Section V-C wasused, which now involves G φ ( k , (cid:126)r, (cid:126)r (cid:48) ) and G A ( k , (cid:126)r, (cid:126)r (cid:48) ) askernels due to the layered background medium. A delta-gapport [29] was employed to excite the structure, with the portlocation marked in Fig. 13, which also shows the computedelectric surface current density at GHz. Frequencies from to GHz were considered, over which the structure spanselectrical sizes between . and wavelengths. The magnitudeand phase of the scattering parameter S are shown inthe top and bottom panels of Fig. 14, respectively. In bothcases, the proposed method is in excellent agreement withthe conventional AIM over the entire frequency range. Thisdemonstrates that the proposed method is also capable ofaccurately handling electrically large problems where the MGFis used to model layered background media. The top panel ofFig. 15 again demonstrates the relative efficacy of the con-straint preconditioner involving only frequency-independentself-term entries, while the bottom panel of Fig. 15 shows theper-frequency CPU time improvement of the proposed methodover the conventional AIM. Table II reports CPU times for thetwo methods, and shows that the near-region direct integrationstep in the proposed method is significantly faster than in theconventional AIM. A total speed-up of . × was obtained withthe proposed method. In this case, the iterative solve time isa significant contributor to the total simulation time at eachfrequency, and is the limiting factor in the obtained speed-up. E. Stacked Inductor Array
As a final example, we consider an on-chip stacked inductorarray [3] consisting of × × copper coils, shown inFig. 16. Each individual element is a × scaled version ofthe geometry described in [58]. The structure is embedded ina six-layer substrate with the configuration shown in Table I.The first inductor layer starts at z = 45 µ m, and the second at z = 68 µ m, where z = 0 is taken at the interface between thelowermost layer and the PEC backing. A mesh with , triangles and , edges was generated, and an AIM gridwith × × points was used with n = 2 . The AEFIEformulation with the same constraint preconditioner as in theprevious two examples was used, and the presence of skineffect was modeled with the surface impedance boundary EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 12
Fig. 13: Geometry and electric surface current density at GHz for the dipole array in Section V-D. . . . . . . . . f (GHz) − − | S | ( d B ) Conventional AIMProposed . . . . . . . . f (GHz) − S ( d e g r ee s ) Conventional AIMProposed
Fig. 14: Accuracy validation for the dipole array in Sec-tion V-D. Top panel: magnitude of S . Bottom panel: phaseof S . . . . . . . . . f (GHz)0200400 I t e r a t i o n s Conventional AIMProposed . . . . . . . . f (GHz)0200400600 C P U t i m e ( s ) Conventional AIMProposed
Fig. 15: Performance comparison for the dipole array inSection V-D. Top panel: iterations required for convergence.Bottom panel: CPU time per frequency.condition (SIBC) [59]. The MGF involving G φ ( k , (cid:126)r, (cid:126)r (cid:48) ) and G A ( k , (cid:126)r, (cid:126)r (cid:48) ) was used to model the background medium. A Fig. 16: Geometry and electric surface current density at GHz for the inductor array in Section V-E. f (GHz) − − | S | ( d B ) | S || S || S || S || S || S | Conventional AIMProposed f (GHz) − S ( d B ) S S S S S S Conventional AIMProposed
Fig. 17: Accuracy validation for the dipole array in Sec-tion V-D. Top panel: magnitude of S . Bottom panel: phaseof S . f (GHz)50100150 I t e r a t i o n s Conventional AIMProposed f (GHz)10 T o t a l t i m e ( s ) Conventional AIMProposed
Fig. 18: Performance comparison for the inductor array inSection V-E. Top panel: iterations required for convergence.Bottom panel: CPU time per frequency.lumped port excitation model was used [22], [60], with portdefinitions shown in Fig. 16, along with the computed electricsurface current density at GHz.The top and bottom panels of Fig. 17 show, respectively, themagnitude and phase of S parameters corresponding to reflec- EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 13 tion, transmission, and coupling, for the proposed method andthe conventional AIM. The proposed method is in excellentagreement in both magnitude and phase over the entire fre-quency range considered. Furthermore, the top panel of Fig. 18confirms the effectiveness of the frequency-independent con-straint preconditioner for all frequencies, while the bottompanel of Fig. 18 demonstrates the significant improvement inCPU time per frequency. The near-region direct integrationtimes are shown in Table II, and the proposed method leadsto a speed-up of . × in total, compared to the conventionalAIM. VI. D ISCUSSION
Through various numerical examples, we demonstrated thatthe proposed method greatly speeds up frequency sweeps inSIE-based electromagnetic solvers without loss of accuracy.Specifically, we emphasize the following advantages arisingfrom the use of the proposed technique: • The cost-per-frequency of the matrix-fill and precorrec-tion steps is eliminated. Instead, these steps need to beperformed only once at the beginning of the frequencysweep. • The computation of near-region interactions is greatlysimplified, since they only involve static kernels. Forlayered media, the implication is that expensive Som-merfeld integration or approximation techniques suchas the DCIM are avoided when computing near-regioninteractions. • For complex geometries, there is usually a tradeoff be-tween mesh quality and CPU time, since a high qualitymesh may require more elements, which causes a sig-nificant increase in the CPU time associated with near-region interactions. The cost of using high-quality meshesis significantly alleviated with the proposed method, sincenear-region interactions are computed and precorrectedonly once. • One may precompute and store to disk the static near-region matrix for a given structure. The matrix could beloaded and reused whenever needed, for example whenfurther refinement is needed near a resonance, or whenthe frequency range must be extended.VII. C
ONCLUSION
A novel technique is proposed for the accelerated electro-magnetic modeling of large problems in homogeneous andlayered media, based on an extended adaptive integral method(AIMx). The kernels of the SIE operators are decomposedinto frequency-independent and frequency-dependent parts,such that the singularity of the kernel is contained entirelyin the frequency-independent part. The frequency-dependentpart is a smooth function which can be accurately interpo-lated, as shown with analytical and numerical error analyses.Therefore, matrix-vector products involving the frequency-dependent terms are fully accelerated with the AIM in boththe near- and far-region. Only the near-region computationsassociated with the frequency-independent term need to becomputed accurately with direct integration, which can be performed in advance of a frequency sweep. This leads tosignificantly faster frequency sweeps. The efficiency of theproposed AIMx is demonstrated through diverse numericalexamples, leading to total speed-ups ranging from × to × .R EFERENCES[1] J.-M. Jin,
The Finite Element Method in Electromagnetics , 3rd ed.Hoboken, NJ, USA: Wiley, 2014.[2] A. E. Ruehli, G. Antonini, and L. Jiang,
Circuit Oriented Electromag-netic Modeling Using the PEEC Techniques . IEEE Press, 2017.[3] T. Moselhy, X. Hu, and L. Daniel, “pFFT in FastMaxwell: A fastimpedance extraction solver for 3D conductor structures over substrate,”in
Proc. Conf. Des., Automat. Test , 2007.[4] S. Omar and D. Jiao, “A new volume integral formulation for broadband3-D circuit extraction in inhomogeneous materials with and withoutexternal electromagnetic fields,”
IEEE Trans. Microw. Theory Tech. ,vol. 61, no. 12, pp. 4302–4312, Dec. 2013.[5] A. Poggio and E. Miller, “Integral equation solutions of three-dimensional scattering problems,” in
Computer Techniques for Elec-tromagnetics , ser. International Series of Monographs in ElectricalEngineering. Pergamon, 1973, pp. 159 – 264.[6] W. Chai and D. Jiao, “Direct matrix solution of linear complexity forsurface integral-equation-based impedance extraction of complicated 3-D structures,”
Proc. IEEE , vol. 101, no. 2, pp. 372–388, Jun. 2013.[7] Z.-G. Qian and W. C. Chew, “Fast full-wave surface integral equa-tion solver for multiscale structure modeling,”
IEEE Trans. AntennasPropag. , vol. 57, no. 11, pp. 3594–3601, Nov. 2009.[8] T. Xia, H. Gan, M. Wei, W. C. Chew, H. Braunisch, Z. Qian, K. Ayg¨un,and A. Aydiner, “An enhanced augmented electric-field integral equationformulation for dielectric objects,”
IEEE Trans. Antennas Propag. ,vol. 64, no. 6, pp. 2339–2347, Jun. 2016.[9] M. Huynen, K. Y. Kapusuz, X. Sun, G. Van der Plas, E. Beyne,D. De Zutter, and D. Vande Ginste, “Entire domain basis functionexpansion of the differential surface admittance for efficient broadbandcharacterization of lossy interconnects,”
IEEE Trans. Microw. TheoryTech. , vol. 68, no. 4, pp. 1217–1233, Jan. 2020.[10] S. Sharma and P. Triverio, “SLIM: A well-conditioned single-sourceboundary element method for modeling lossy conductors in layeredmedia,”
IEEE Antennas Wireless Propag. Lett. , 2020 (submitted, axXiv:2007.07378).[11] L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,”
Journal of Computational Physics , vol. 73, no. 2, pp. 325 – 348, 1987.[12] N. Engheta, W. D. Murphy, V. Rokhlin, and M. S. Vassiliou, “Thefast multipole method (FMM) for electromagnetic scattering problems,”
IEEE Trans. Antennas Propag. , vol. 40, no. 6, pp. 634–641, Jun. 1992.[13] J. Song, C.-C. Lu, and W. C. Chew, “Multilevel fast multipole algorithmfor electromagnetic scattering by large complex objects,”
IEEE Trans.Antennas Propag. , vol. 45, no. 10, pp. 1488–1493, Oct. 1997.[14] L. J. Jiang and W. C. Chew, “A mixed-form fast multipole algorithm,”
IEEE Trans. Antennas Propag. , vol. 53, no. 12, pp. 4145–4156, Dec.2005.[15] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, “AIM: Adaptiveintegral method for solving large-scale electromagnetic scattering andradiation problems,”
Radio Sci. , vol. 31, no. 5, pp. 1225–1251, Sep.1996.[16] K. A. Michalski and J. R. Mosig, “Multilayered media Green’s functionsin integral equation formulations,”
IEEE Trans. Antennas Propag. ,vol. 45, no. 3, pp. 508–519, Mar. 1997.[17] F. Ling, J. Song, and J.-M. Jin, “Multilevel fast multipole algorithmfor analysis of large-scale microstrip structures,”
IEEE Microw. GuidedWave Lett. , vol. 9, no. 12, pp. 508–510, Dec. 1999.[18] B. Hu and W. C. Chew, “Fast inhomogeneous plane wave algorithm forscattering from objects above the multilayered medium,”
IEEE Trans.Geosci. Remote Sens. , vol. 39, no. 5, pp. 1028–1038, May 2001.[19] V. Okhmatovski, M. Yuan, I. Jeffrey, and R. Phelps, “A three-dimensional precorrected FFT algorithm for fast method of momentssolutions of the mixed-potential integral equation in layered media,”
IEEE Trans. Microw. Theory Tech. , vol. 57, no. 12, pp. 3505–3517,Dec. 2009.[20] F. Ling, C.-F. Wang, and J.-M. Jin, “An efficient algorithm for ana-lyzing large-scale microstrip structures using adaptive integral methodcombined with discrete complex-image method,”
IEEE Trans. Microw.Theory Tech. , vol. 48, no. 5, pp. 832–839, May 2000.
EEE TRANSACTIONS ON ANTENNAS AND PROPAGATION 14 [21] K. Yang and A. E. Yilmaz, “A three-dimensional adaptive integralmethod for scattering from structures embedded in layered media,”
IEEETrans. Geosci. Remote Sens. , vol. 50, no. 4, pp. 1130–1139, Apr. 2012.[22] S. Sharma, U. R. Patel, S. V. Hum, and P. Triverio, “A complete surfaceintegral method for broadband modeling of 3D interconnects in stratifiedmedia,” arXiv e-prints , p. arXiv: 1810.04030, Oct. 2018.[23] Z. Zhu, B. Song, and J. K. White, “Algorithms in FastImp: a fastand wide-band impedance extraction program for complicated 3-Dgeometries,”
IEEE Trans. Comput.-Aided Design Integr. Circuits Syst. ,vol. 24, no. 7, pp. 981–998, Jul. 2005.[24] S. Q. Li, Y. Yu, C. H. Chan, K. F. Chan, and L. Tsang, “A sparse-matrix/canonical grid method for analyzing densely packed intercon-nects,”
IEEE Trans. Microw. Theory Tech. , vol. 49, no. 7, pp. 1221–1228,Jul. 2001.[25] B. J. Fasenfest, F. Capolino, D. R. Wilton, D. R. Jackson, and N. J.Champagne, “A fast MoM solution for large arrays: Green’s functioninterpolation with FFT,”
IEEE Antennas Wireless Propag. Lett. , vol. 3,pp. 161–164, 2004.[26] S. M. Seo and J.-F. Lee, “A fast IE-FFT algorithm for solving PECscattering problems,”
IEEE Trans. Magn. , vol. 41, no. 5, pp. 1476–1479,May 2005.[27] H. G. Wang and C. H. Chan, “The implementation of multilevel Green’sfunction interpolation method for full-wave electromagnetic problems,”
IEEE Trans. Antennas Propag. , vol. 55, no. 5, pp. 1348–1358, May2007.[28] D. Schobert and T. Eibert, “Direct field and mixed potential integralequation solutions by fast Fourier transform accelerated multilevelGreen’s function interpolation for conducting and impedance bound-ary objects,”
Applied Computational Electromagnetics Society Journal ,vol. 26, pp. 1016–1023, Jan. 2011.[29] W. C. Gibson,
The Method of Moments in Electromagnetics . BocaRaton, FL, USA: CRC press, 2014.[30] B. Fasenfest, J. D. Rockway, N. J. Champagne, and R. M. Sharpe,“A generalized fast frequency sweep algorithm for coupled circuit-EMsimulations,” in
IEEE Int. Symp. Antennas Propag. , vol. 4, Sep. 2004,pp. 3944–3947.[31] T. Dhaene, J. Ureel, N. Fache, and D. De Zutter, “Adaptive frequencysampling algorithm for fast and accurate S-parameter modeling ofgeneral planar structures,” in
Proc. of 1995 IEEE MTT-S Int. Microw.Symp. , 1995, pp. 1427–1430 vol.3.[32] R. Lehmensiek and P. Meyer, “Creating accurate multivariate rationalinterpolation models of microwave circuits by using efficient adaptivesampling to minimize the number of computational electromagneticanalyses,”
IEEE Trans. Microw. Theory Tech. , vol. 49, no. 8, pp. 1419–1430, Aug. 2001.[33] G. Antonini, D. Deschrijver, and T. Dhaene, “Broadband rational macro-modeling based on the adaptive frequency sampling algorithm and thepartial element equivalent circuit method,”
IEEE Trans. Electromagn.Compat. , vol. 50, no. 1, pp. 128–137, Feb. 2008.[34] G. Antonini and D. Romano, “An accurate interpolation strategy for fastfrequency sweep of partial element equivalent circuit models,”
IEEETrans. Electromagn. Compat. , vol. 56, no. 3, pp. 653–658, Nov. 2014.[35] D. Colton and R. Kress,
Integral Equation Methods in Scattering Theory .Hoboken, NJ, USA: Wiley, 1983.[36] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering bysurfaces of arbitrary shape,”
IEEE Trans. Antennas Propag. , vol. 30,no. 3, pp. 409–418, May 1982.[37] M. M. Jia, S. Sun, Y. Li, Z. Qian, and W. C. Chew, “Accelerationof perturbation-based electric field integral equations using fast Fouriertransform,”
IEEE Trans. Antennas Propag. , vol. 64, no. 10, pp. 4559–4564, Oct. 2016.[38] L. Lombardi, G. Antonini, and A. E. Ruehli, “Analytical evaluation ofpartial elements using a retarded Taylor series expansion of the Green’sfunction,”
IEEE Trans. Microw. Theory Tech. , vol. 66, no. 5, pp. 2116–2127, May 2018.[39] K. Konno, Q. Chen, and R. J. Burkholder, “Fast computation of layeredmedia Green’s function via recursive Taylor expansion,”
IEEE AntennasWireless Propag. Lett. , vol. 16, pp. 1048–1051, Oct. 2017.[40] S. Sharma, U. R. Patel, P. Triverio, “Accelerated electromagnetic analysisof interconnects in layered media using a near-field series expansion ofthe Greens function,” in
IEEE Conf. Electr. Perform. Electron. Packag.and Syst. , San Jose, CA, Oct. 2018.[41] W.-B. Ewe, L.-W. Li, Q. Wu, and M.-S. Leong, “Preconditionersfor adaptive integral method implementation,”
IEEE Trans. AntennasPropag. , vol. 53, no. 7, pp. 2346–2350, Jul. 2005. [42] O. Wiedenmann and T. F. Eibert, “The effect of near-zone precondi-tioning on electromagnetic integral equations of first and second kind,”
Advances in Radio Science , vol. 11, pp. 61–65, 2013.[43] W. C. Chew,
Waves and Fields in Inhomogeneous Media . Hoboken,NJ, USA: Wiley, 1999.[44] K. A. Michalski and D. Zheng, “Electromagnetic scattering and radiationby surfaces of arbitrary shape in layered media. i. theory,”
IEEE Trans.Antennas Propag. , vol. 38, no. 3, pp. 335–344, Mar. 1990.[45] K. A. Michalski and J. R. Mosig, “Efficient computation of Sommerfeldintegral tails - methods and algorithms,”
Journal of ElectromagneticWaves and Applications , vol. 30, no. 3, pp. 281–317, 2016.[46] M. I. Aksun, “A robust approach for the derivation of closed-formgreen’s functions,”
IEEE Trans. Microw. Theory Tech. , vol. 44, no. 5,pp. 651–658, May 1996.[47] E. Simsek, Q. H. Liu, and B. Wei, “Singularity subtraction for evaluationof green’s functions for multilayer media,”
IEEE Trans. Microw. TheoryTech. , vol. 54, no. 1, pp. 216–225, Jan. 2006.[48] K. Yang and A. E. Yilmaz, “Comparison of precorrected FFT/adaptiveintegral method matching schemes,”
Microwave and Optical TechnologyLetters , vol. 53, no. 6, pp. 1368–1372, Jun. 2011.[49] M. Abramowitz and I. A. Stegun,
Handbook of Mathematical Functionswith Formulas, Graphs, and Mathematical Tables , 10th ed. New York:Dover, 1964.[50] S. Sharma and P. Triverio, “A fully-accelerated surface integral equationmethod for the electromagnetic modeling of arbitrary objects,”
IEEETrans. Antennas Propag. , 2020 (submitted, arXiv: 2003.11679).[51] The Sage Developers,
SageMath, the Sage Mathematics Software System(Version 9.1) , 2020, .[52] S. Balay, S. Abhyankar M. F. Adams, J. Brown, P. Brune, K. Buschel-man, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley,D. A. May, L. C. McInnes R. T. Mills, T. Munson, K. Rupp, P. SananB. F. Smith, S. Zampini, H. Zhang.
PETSc Web Page
SIAM J. ScientificStatistical Comput. , vol. 7, no. 3, pp. 856–869, 1986.[54] A. C. Woo, H. T. G. Wang, M. J. Schuh, and M. L. Sanders, “EMprogrammer’s notebook-benchmark radar targets for the validation ofcomputational electromagnetics programs,”
IEEE Antennas Propag.Mag. , vol. 35, no. 1, pp. 84–89, Feb. 1993.[55] U. R. Patel, P. Triverio, and S. V. Hum, “A macromodeling approach toefficiently compute scattering from large arrays of complex scatterers,”
IEEE Trans. Antennas Propag. , vol. 66, no. 11, pp. 6158–6169, Nov.2018.[56] C. Geaney, M. Hosseini, and S. Hum, “A cascaded polarizerreflectar-ray for independent dual circular polarization control,”
ESA AntennaWorkshop on Antennas and RF Systems for Space Science, Noordwijk,Netherlands, Oct. 2017 , Oct. 2017.[57] B. G. Duffley, G. A. Morin, M. Mikavica, and Y. M. M. Antar, “A wide-band printed double-sided dipole array,”
IEEE Trans. Antennas Propag. ,vol. 52, no. 2, pp. 628–631, Feb. 2004.[58] U. R. Patel, S. Sharma, S. Yang, S. V. Hum, and P. Triverio, “Full-waveelectromagnetic characterization of 3D interconnects using a surfaceintegral formulation,” in
IEEE Conf. Electr. Perform. Electron. Packag.and Syst. , San Jose, CA, Oct. 2017.[59] S. V. Yuferev and N. Ida,
Surface Impedance Boundary Conditions: AComprehensive Approach . Boca Raton, FL, USA: CRC Press, 2009.[60] Y. Wang, D. Gope, V. Jandhyala, and C.-J. R. Shi, “GeneralizedKirchoff’s current and voltage law formulation for coupled circuit-electromagnetic simulation with surface integral equations,”