A Chebyshev-based High-order-accurate Integral Equation Solver for Maxwell's Equations
11 A Chebyshev-based High-order-accurate IntegralEquation Solver for Maxwells Equations
Jin Hu,
Student Member, IEEE and Constantine Sideris,
Member, IEEE
Abstract —This paper introduces a new method for discretizingand solving integral equation formulations of Maxwell’s equa-tions which achieves spectral accuracy for smooth surfaces. Theapproach is based on a hybrid Nystr¨om-collocation method usingChebyshev polynomials to expand the unknown current densitiesover curvilinear quadrilateral surface patches. As an example,the proposed strategy is applied to Magnetic Field IntegralEquation (MFIE) and the N-M¨uller formulation for scatteringfrom metallic and dielectric objects respectively. The convergenceis studied for several different geometries, including spheres,cubes, and complex NURBS geometries imported from CADsoftware, and the results are compared against a commercialMethod-of-Moments solver using RWG basis functions.
Index Terms —Integral equations, high-order accuracy, N-M¨uller formulation, spectral methods, scattering.
I. I
NTRODUCTION D UE to the lack of analytical solutions for anything but thesimplest problems [1], efficient and accurate numericalmethods for solving Maxwell’s equations are crucial for aplethora of engineering applications today, including antennas,microwave devices, and nanophotonic structures. A recentresurgence in inverse design approaches [2], which involve theautomated design of novel electromagnetic structures givena set of desired performance metrics and design constraintsand require accurate field and gradient information at eachiteration, further highlights the need for fast Maxwell solvers.Although finite difference [3] and finite element methods [4]are popular approaches due to their relative ease of implemen-tation, they suffer from several major drawbacks: poor con-vergence due to finite difference approximations or low-orderbasis functions, significant numerical dispersion due to relyingon local discrete differentiation, and they are often impracticalfor large problems due to their volumetric nature. On the otherhand, boundary integral equation (BIE) formulations have beenshown to be very effective in situations containing scattererswith small surface area to volume ratios due to only solvingfor unknowns on surfaces rather than volumes. Recently,BIEs have been successfully applied towards the modelingand optimization of nanophotonic devices in two dimensions,showing significant improvements in speed and accuracy overfinite difference based methods [5]. The majority of presentday implementations of BIE methods rely on discretization ofobjects via flat triangular discretizations, expanding unknownsurface current densities with Rao-Wilton-Glisson (RWG)
This authors gratefully acknowledge support by the National Science Foun-dation (1849965) and the Air Force Office of Scientific Research (FA9550-20-1-0087).J. Hu and C. Sideris are with the Department of Electrical and ComputerEngineering, University of Southern California, Los Angeles, CA 90089, USA(e-mails: [email protected], [email protected]). basis functions, and utilizing the Method of Moments (MoM)for solving the resulting integral formulation [6]. Unfortu-nately, RWG functions are first order and cannot accuratelyapproximate complex surface current distributions withoutvery fine meshing. This often leads to poor convergence andconditioning of these implementations and significant researcheffort has gone into selecting alternative basis functions fortesting or expansion in an attempt to improve performance [7].[8] presented an alternative approach which achieves high-order accuracy by utilizing a Nystr¨om method and discretizingthe integrals on the basis of local coordinate charts togetherwith fixed and floating partitions of unity. While effective, theapproach of [8] relies on overlapping parameterized patcheswhich can both increase the number of unknowns as wellas significantly complicate the generation of surface meshes.Recently, [9] demonstrated a new high-order solution strategyfor acoustic scattering problems based on non-overlappingparametric curvilinear patches. The method presented in [9]discretizes the unknowns on each patch on a Chebyshev grid,approximating the unknown surface densities using Chebyshevpolynomials. A spectrally accurate Fej`er quadrature rule isused for evaluating far interactions, and a Cartesian changeof variables is used to cancel the singularity of the integralsassociated with local and near interactions, leading to high-order accuracy in the numerical evaluation of both the singularand near-singular integrals.In this work, we extend the methods presented in [9] forthe scalar Helmholtz equation to the numerical solution ofthe fully-vectorial Maxwell case. In order to demonstrate thegenerality of the approach, we consider scattering from bothPerfect Electrical Conductor (PEC) and dielectric objects.We focus on the solution of the MFIE formulation [10]for metallic objects and the N-M¨uller formulation [11] fordielectric objects due to their superior conditioning properties,although we remark that all of the methods presented inthis work can readily be extended to the Electric and Com-bined Field Integral Equations (EFIE/CFIE) and other integralequation formulations designed for dielectric objects, suchas the Poggio-Chang-Miller-Harrington-Wu-Tsai (PCMHWT)formulation [12]. This paper is organized as follows. InSection II, we briefly review the MFIE and the N-M¨ullerformulations. In Section III, we present the proposed high-order-accurate Chebyshev-based Boundary Integral Equation(CBIE) approach for discretizing the integral formulations.Finally, numerical results are presented in Section IV whichevaluate the performance of the CBIE method by compar-ing the numerical solutions of plane-wave scattering from aPEC/dielectric sphere against analytical Mie-series solutions,as well as solving a PEC/dielectric cube for which no closed- a r X i v : . [ phy s i c s . c o m p - ph ] J u l form solutions exist. The accuracy is also compared againsta commercial RWG-based MoM solver. Finally, we presentresults for scattering from two complex NURBS parametrizedgeometries generated by commerical CAD software.II. I NTEGRAL E QUATION F ORMULATIONS
A. Magnetic Field Integral Equation Formulation for ClosedMetallic Scatterers
We consider the problem of computing the scattered electricand magnetic fields ( E scat , H scat ) that result due to an incidentfield excitation (cid:0) E inc , H inc (cid:1) impinging on the surface Γ of aclosed perfect metallic object D as illustrated in Fig. 1(a).Based on the Stratton-Chu formulas [13], Electric and Mag-netic Field Integral Equations (EFIE/MFIE) can be derivedwhich express the scattered electric and magnetic fields interms of the physical current J = ˆn × H on the surface of aperfect metallic conducting object. Although either the EFIE,the MFIE, or a linear combination of the two can be usedto solve for the scattered fields due to an incident excitation,only the MFIE is considered in this work due to its superiorconditioning properties as a result of the nature of Fredholmintegral equations of the second kind [13]. The classical MFIEcan be expressed as, J K J = ˆn × H inc (1)where K is the operator: K [ a ] ( r ) = ˆn ( r ) × (cid:90) Γ a ( r (cid:48) ) × ∇ G ( r − r (cid:48) ) dσ ( r (cid:48) ) (2)Note that ∇ denotes the gradient with respect to the coor-dinates of observation points r , G corresponds to the freespace scalar Green’s function: G ( r − r (cid:48) ) = exp ( − jk | r − r (cid:48) | ) π | r − r (cid:48) | with wavenumber k = 2 π/λ , and ˆn denotes the outwardlypointing surface normal. B. N-M¨uller Formulation for Dielectric Scatterers
The second scenario that we consider is scattering froma penetrable dielectric object D with a permittivity (cid:15) d anda permeability µ d embedded in a homogeneous backgroundmedium characterized by permittivity (cid:15) e and permeability µ e in the presence of an incident field excitation (cid:0) E inc , H inc (cid:1) .As shown in Fig. 1(b), since the object is now penetrable,the incident fields lead to scattered fields outside the object, ( E scat , H scat ) , as well as transmitted fields inside, ( E t , H t ) .Equivalent electric and magnetic current densities can thenbe defined based on the boundary tangential magnetic andelectric fields respectively across the dielectric interface as: J = ˆn × ( H inc + H scat ) = ˆn × H t and M = ( E inc + E scat ) × ˆn = E t × ˆn on the surface Γ of D . By invoking the Stratton-Chuformula for the electric and magnetic fields outside of theobject and crossing with the normal vector ˆn , we obtain: M K e M − η e T e J = − ˆn × E inc (3) J K e J + 1 η e T e M = ˆn × H inc (4) (a)(b) Fig. 1: (a) EM scattering from a closed PEC object. (b) EMscattering from a closed penetrable dielectric object.where the K e and T e operators are defined as: K e [ a ] ( r ) = ˆn ( r ) × (cid:90) Γ a ( r (cid:48) ) × ∇ G e ( r − r (cid:48) ) dσ ( r (cid:48) ) (5) T e [ a ] ( r ) = T se [ a ] ( r ) + T he [ a ] ( r ) (6) T se [ a ] ( r ) = jk e ˆn ( r ) × (cid:90) Γ a ( r (cid:48) ) G e ( r − r (cid:48) ) dσ ( r (cid:48) ) (7) T he [ a ] ( r ) = jk e ˆn ( r ) × (cid:90) Γ ∇ G e ( r − r (cid:48) ) ∇ (cid:48) s · a ( r (cid:48) ) dσ ( r (cid:48) ) (8)where the subscript “e” in the operators indicates the exteriormedium, which has wavenumber: k e = 2 π/λ e and impedance: η e = (cid:112) µ e /(cid:15) e .Similarly, another set of integral equations can be obtainedfor the transmitted fields ( E t , H t ) inside the object: M − K d M + η d T d J = (9) J − K d J − η d T d M = (10)where the K d and T d operators are defined in the same manneras K e and K d , except the “d” denotes the interior mediumwith corresponding wavenumber: k d = 2 π/λ d and impedance: η d = (cid:112) µ d /(cid:15) d Fig. 2: The mapping from square [-1,1]x[-1,1] in parameterdomain to a patch on a sphere in Cartesian coordinates.(3), (4), (9), and (10) give four equations for two unknowns ( J , M ) . They can be linearly combined as follows to reducethe system to two independent equations: α (3) + α (9) β (4) + β (10) (11)Choosing α = (cid:15) e , α = (cid:15) d , β = µ e , β = µ d results in theclassical N-M¨uller formulation, which has the advantage that itcompletely cancels the singular terms arising from the gradientof the Green’s function in the T he and T hd operators [14]. Thecombined system in matrix form is thus: (cid:20) (cid:15) e K e − (cid:15) d K d + (cid:15) e + (cid:15) d I − ( MT s + MT h ) MT s + MT h µ e K e − µ d K d + µ e + µ d I (cid:21) (cid:20) MJ (cid:21) = (cid:20) − (cid:15) e ˆn × E inc µ e ˆn × H inc (cid:21) (12)where I is the identity operator, MT s and MT h are definedas MT s [ a ] ( r ) = ( √ µ e (cid:15) e T se − √ µ d (cid:15) d T sd ) [ a ] ( r )= jω ˆn ( r ) × (cid:90) Γ a ( r (cid:48) )( k e G e − k d G d ) dσ ( r (cid:48) ) (13) MT h [ a ] ( r ) = ( √ µ e (cid:15) e T he − √ µ d (cid:15) d T hd ) [ a ] ( r )= jω ˆn ( r ) × (cid:90) Γ ( ∇ G e − ∇ G d ) ∇ (cid:48) s · a ( r (cid:48) ) dσ ( r (cid:48) ) (14)The difference of the hypersingular operators T h , MT h doesnot have any singularity due to exact cancellation of thesingular terms.III. C HEBYSHEV - BASED B OUNDARY I NTEGRAL E QUATION A PPROACH
A. Representation of Geometries/Densities
In order to solve (1) or (12), the surface Γ is first divided intoa number ( M ) of non-overlapping curvilinear patches Γ p , p =1 , , ..., M . We opt to use parametric curvilinear patches,rather than a more common flat triangular discretization, due totheir ability to accurately represent curved surfaces with coarsediscretizations. For each of these patches, a U V mapping isused to map from the standard [ − , × [1 , square in U V space to the corresponding parameterized surface in Cartesiancoordinates as illustrated in Fig. 2. Defining the position vectoron Γ p as r = r p ( u, v ) = ( x p ( u, v ) , y p ( u, v ) , z p ( u, v )) , we can define the tangential covariant basis vectors and surface normalon Γ p as a pu = ∂ r p ( u, v ) ∂u , a pv = ∂ r p ( u, v ) ∂v , ˆn p = a pu × a pv || a pu × a pv || (15)Thus, the vector triplet ( a pu , a pv , ˆn p ) forms a local conformalreference frame at each point on Γ p . The metric tensor isdefined as G p = (cid:20) g puu g puv g pvu g pvv (cid:21) (16)where g pij = a pi · a pj and thus we have surface element ds = (cid:112) | G p | dudv on Γ p where | G p | is the determinant of G p . Wecan now represent the surface current density on Γ p as J ( r p ( u, v )) = J p ( u, v )= 1 (cid:112) | G p ( u, v ) | ( J p,u ( u, v ) a pu ( u, v ) + J p,v ( u, v ) a pv ( u, v )) (17) M ( r p ( u, v )) = M p ( u, v )= 1 (cid:112) | G p ( u, v ) | ( M p,u ( u, v ) a pu ( u, v ) + M p,v ( u, v ) a pv ( u, v )) (18)for p = 1 , ..., M , where J p,u (resp. M p,u ) and J p,v (resp. M p,v ) are scalar functions representing the contravariant com-ponents of the surface current density J (resp. M ) on the p th patch normalized by the metric tensor, (cid:112) | G p | . The densitiesare normalized by the surface element in order to simplifythe numerical computation of their divergence (see [13, sec.6.2.5]). Due to their desirable spectral convergence propertiesfor approximating smooth functions, we utilize Chebyshevpolynomials to discretize the surface current densities: J p,a = N pv − (cid:88) m =0 N pu − (cid:88) n =0 γ p,an,m T n ( u ) T m ( v ) , for a = u, v (19) M p,a = N pv − (cid:88) m =0 N pu − (cid:88) n =0 ζ p,an,m T n ( u ) T m ( v ) , for a = u, v (20)where the Chebyshev coefficients γ p,jn,m and ζ p,jn,m can becomputed from the values of the densities on Chebyshevnodes, γ p,an,m = α n α m N pu N pv N pv − (cid:88) k =0 N pu − (cid:88) l =0 J p,a ( u l , v k ) T n ( u l ) T m ( v k ) , (21) ζ p,an,m = α n α m N pu N pv N pv − (cid:88) k =0 N pu − (cid:88) l =0 M p,a ( u l , v k ) T n ( u l ) T m ( v k ) , (22)based on the discrete-orthogonality property of Chebyshevpolynomials [15], where α n = (cid:40) , n = 02 , n (cid:54) = 0 (23)Therefore, only the unknowns at the Chebyshev nodes (37)are required to represent the continuous scalar densities J p,a and M p,a over the whole patch Γ p , where a can be either u or v . In our specific implementation, these unknowns are orderedin vector form as: J p = J p,u ( u , v ) ... J p,u ( u N pu − , v N pv − ) J p,v ( u , v ) ... J p,v ( u N pu − , v N pv − ) , (24) M p = M p,u ( u , v ) ... M p,u ( u N pu − , v N pv − ) M p,v ( u , v ) ... M p,v ( u N pu − , v N pv − ) , (25)Here, ( u l , v k ) , l = 0 , , ..., N pu − , k = 0 , , ..., N pv − arethe Cartesian product of the Chebyshev nodes in the u and v direction as given in (37) for a discretization using N pu · N pv points on Γ p . B. Discretization of operators
We now turn our attention towards discretization of the K / K e / K d , MT s and MT h operators. We will begin bydiscretizing the K operator first. Clearly, any integral over Γ can be split into the sum of integrals over each of the M patches, K [ J ] ( r ) = M (cid:88) p =1 K [ J p ] K [ J p ] ( r ) = ˆn ( r ) × (cid:90) Γ p J p ( r (cid:48) ) × ∇ G ( r − r (cid:48) ) dσ ( r (cid:48) )= ˆn ( r ) × (cid:90) − (cid:90) − ( J p,u ( u, v ) a pu ( u, v ) + J p,v ( u, v ) a pv ( u, v )) × ∇ G ( r − r p ( u, v )) dudv (26)Note that the (cid:112) | G p ( u, v ) | in the denominator of the expansion(17) for J cancels with the Jacobian (cid:112) | G p ( u, v ) | that appearsin the integral. In its current form, (26) contains the hyper-singular kernel ∇ G ; however, it can be manipulated using theBAC-CAB vector identity into K [ J p ] ( r ) = (cid:90) − (cid:90) − J p,u ( u, v )( a pu ( u, v ) ∂G ( r − r p ( u, v )) ∂ ˆn ( r ) − ∇ G ( r − r p ( u, v )) ˆn ( r ) · a pu ( u, v )) dudv + (cid:90) − (cid:90) − J p,v ( u, v )( a pv ( u, v ) ∂G ( r − r p ( u, v )) ∂ ˆn ( r ) − ∇ G ( r − r p ( u, v )) ˆn ( r ) · a pv ( u, v )) dudv (27)which is weakly singular since ˆn ( r ) · a pu,v approaches 0 as r p ( u, v ) → r . Substituting (27) into (1), we must obtain N = 2 (cid:80) Mp =1 N pu N pv linearly independent equations in orderto obtain a uniquely solvable linear system for approximating J on Γ . This is achieved by using a collocation method andtesting (1) at same points as the unknowns. Since J is a vector function with two unknown contravariant components, wemust test the MFIE at each point with two linearly independentvectors. The natural choice for this is the set of normalizedcovariant basis vectors √ G p a p,u and √ G p a p,v where thecovariant basis vectors a p,u and a p,v are defined via theorthogonality relation a p,a · a pb = (cid:40) a = b a (cid:54) = b (28)since dotting (1) with them results in expanding the right handside ˆn × H inc using the basis vectors as the unknown densities.We can now define the linear system: I + K K . . . K M K I + K . . . K M ... ... . . . ... K M K M . . . I + K MM J J ... J M = H inc H inc ... H M inc (29)where H p inc = − a pv · H p, inc ( u , v ) ... − a pv · H p, inc ( u N pu − , v N pv − ) a pu · H p, inc ( u , v ) ... a pu · H p, inc ( u N pu − , v N pv − ) (30)represents the incident magnetic field on the p th patch and J p , p = 1 , , .., M is given in (24). Note the vector identitiesused in arriving at the expression in (30) are: √ G p a p,u · ˆn p × H p, inc = a p,u · (( a pu × a pv ) × H p, inc )= H p, inc · ( a p,u × ( a pu × a pv ))= H p, inc · ( a pu ( a p,u · a pv ) − a pv ( a p,u · a pu )) = − a pv · H p, inc (31)Similarly, √ G p a p,v · ˆn p × H p, inc = a pu · H p, inc . Matrix block K qp represents contributions of the appropriately discretized K operator from the densities of the patch p to the target pointson patch q and consists of the individual sub-blocks: K qp = (cid:18) K qpuu K qpuv K qpvu K qpvv (cid:19) (32)For the operators used in the N-M¨uller formulation, thematrix blocks corresponding to the K e and K d operator canbe obtained in exactly the same way as K operator by simplyreplacing the wavenumber k in the Green’s function in (27)with k e and k d respectively. The integral of the MT s and MT h operators can also be split over each patch in as similarway as the K operator: MT s [ J ] ( r ) = M (cid:88) p =1 MT s [ J p ] MT s [ J p ] ( r ) = jω ˆn ( r ) × (cid:90) Γ p J p ( r (cid:48) )( k e G e − k d G d ) dσ ( r (cid:48) )= jω ˆn ( r ) × (cid:90) − (cid:90) − ( J p,u ( u, v ) a pu ( u, v ) + J p,v ( u, v ) a pv ( u, v ))( k e G e ( r − r p ( u, v )) − k d G d ( r − r p ( u, v ))) dudv (33) MT h [ J ] ( r ) = M (cid:88) p =1 MT h [ J p ] MT h [ J p ] ( r ) = jω ˆn ( r ) × (cid:90) Γ p ( ∇ G e − ∇ G d ) ∇ (cid:48) s · J p ( r (cid:48) ) dσ ( r (cid:48) )= jω ˆn ( r ) × (cid:90) − (cid:90) − ( ∇ G e − ∇ G d )( ∂J p,u ∂u + ∂J p,v ∂v ) dudv = jω ˆn ( r ) × (cid:90) − (cid:90) − ( ∇ G e ( r − r p ( u, v )) − ∇ G d ( r − r p ( u, v )))( N pv − (cid:88) m =0 N pu − (cid:88) n =0 γ p,un,m T (cid:48) n ( u ) T m ( v ) + N pv − (cid:88) m =0 N pu − (cid:88) n =0 γ p,vn,m T n ( u ) T (cid:48) m ( v )) dudv = N pv − (cid:88) m =0 N pu − (cid:88) n =0 jω ˆn ( r ) × (cid:90) − (cid:90) − ( γ p,un,m T (cid:48) n ( u ) T m ( v )+ γ p,vn,m T n ( u ) T (cid:48) m ( v ))( ∇ G e ( r − r p ( u, v )) − ∇ G d ( r − r p ( u, v ))) dudv (34)The partial derivative of the densities can be readily computedby taking the derivative of the corresponding Chebyshevpolynomials. After the substitution of (33) and (34) into (12)with the expansion defined in (19) and (20), testing (12) at thesame collocation points as the unknowns results in the linearsystem: (cid:20) (cid:15) e K e − (cid:15) d K d + (cid:15) e + (cid:15) d I − ( M T s + M T h ) M T s + M T h µ e K e − µ d K d + µ e + µ d I (cid:21) (cid:20) MJ (cid:21) = (cid:20) − (cid:15) e E inc µ e H inc (cid:21) (35)The block in E inc corresponding to the incident electric fieldon the p th patch can be obtained in the same manner as in(31): E p inc = − a pv · E p, inc ( u , v ) ... − a pv · E p, inc ( u N pu − , v N pv − ) a pu · E p, inc ( u , v ) ... a pu · E p, inc ( u N pu − , v N pv − ) (36)The counterpart H p inc is defined in (30). The matrices K e , K d , M T s and M T h all have the same block structure arrangedby patches as indicated in (29) and (32) for the matrix K . Asuitable numerical integration strategy must now be chosen forevaluating the necessary operators to compute the above matrixsub-blocks. In the following two subsections, we will detailour approach for dealing with the non-adjacent interactions( p (cid:54) = q ) and the singular and near-singular interactions arisingeither when p = q or when p (cid:54) = q , but the target point on q islocated very near to the source patch p . C. Non-adjacent Interactions
The integrals (27), (33) and (34) are smooth for target pointsfar away from the source patch p . Since the current density J / M is discretized on a Chebyshev grid on each patch, we canuse Fej´er’s first quadrature rule to numerically evaluate theseintegrals with high-order accuracy. The quadrature nodes andweights for an order N open rule are given by: x i = cos (cid:18) π i + 12 N (cid:19) , i = 0 , ..., N − (37) w i = 2 N − N/ (cid:88) k =1 k − (cid:18) kπ i + 1 N (cid:19) (38)and the discretized versions of (27), (33) and (34) become(with a = { u, v } and b = { u, v } ): K qpba [ J p,a ] ( u (cid:48) , v (cid:48) ) = N pv − (cid:88) k =0 N pu − (cid:88) l =0 A qpba ( u (cid:48) , v (cid:48) , u l , v k ) (cid:112) | G q ( u (cid:48) , v (cid:48) ) | w l w k J p,a ( u l , v k ) (39) M T s,qpba [ J p,a ] ( u (cid:48) , v (cid:48) ) = N pv − (cid:88) k =0 N pu − (cid:88) l =0 B qpba ( u (cid:48) , v (cid:48) , u l , v k ) (cid:112) | G q ( u (cid:48) , v (cid:48) ) | w l w k J p,a ( u l , v k ) (40) M T h,qpba [ J p,a ] ( u (cid:48) , v (cid:48) ) = N pv − (cid:88) k =0 N pu − (cid:88) l =0 C qpba ( u (cid:48) , v (cid:48) , u l , v k ) (cid:112) | G q ( u (cid:48) , v (cid:48) ) | w l w k ∂J p,a ∂a ( u l , v k ) (41)with A qpba ( u (cid:48) , v (cid:48) , u l , v k ) = a q,b ( u (cid:48) , v (cid:48) ) · a pa ( u l , v k ) ∂G ( r q ( u (cid:48) , v (cid:48) ) − r p ( u l , v k )) ∂ ˆn q ( u (cid:48) , v (cid:48) ) − ˆn q ( u (cid:48) , v (cid:48) ) · a pa ( u l , v k ) a q,b ( u (cid:48) , v (cid:48) ) · ∇ G ( r q ( u (cid:48) , v (cid:48) ) − r p ( u l , v k )) (42) B qpba ( u (cid:48) , v (cid:48) , u l , v k ) = jω a q,b ( u (cid:48) , v (cid:48) ) · ( ˆn q ( u (cid:48) , v (cid:48) ) × a pa ( u l , v k )) (cid:2) k e G e − k d G d (cid:3) ( r q ( u (cid:48) , v (cid:48) ) − r p ( u l , v k )) (43) C qpba ( u (cid:48) , v (cid:48) , u l , v k ) = jω a q,b ( u (cid:48) , v (cid:48) ) · ˆn q ( u (cid:48) , v (cid:48) ) × [ ∇ G e − ∇ G d ] ( r q ( u (cid:48) , v (cid:48) ) − r p ( u l , v k )) (44)where u l and v k are the discretization points on the Chebyshevgrid corresponding to the x i nodes: u l = x l | l = 0 , . . . , N pu − , v k = x k | k = 0 , . . . , N pv − , and w l and w k are thequadrature weights in the u and v directions respectively. D. Singular and Near-singular Interactions
When the observation point ( u (cid:48) , v (cid:48) ) is on the same patch asthe source patch p , the integrals (27), (33) and (34) becomesingular . In order to accurately compute the resulting integralswith high-order accuracy we consider the following smoothingchange of variables [16, Sec. 3.5]: u ( s ) = ξ u (cid:48) ( s ) , v ( t ) = ξ v (cid:48) ( t ) , for − ≤ s, t ≤ (45) Actually, the integral (34) for MT h remains regular due to the M¨ullercancellation and does not require special consideration; however, for simplicitywe treat it in the same way as the other operators in our implementation. where ξ α ( τ ) = (cid:40) α + (cid:16) sgn ( τ ) − απ (cid:17) w ( π | τ | ) , for α (cid:54) = ± α ∓ (cid:0) ± aπ (cid:1) w (cid:0) π (cid:12)(cid:12) τ ∓ (cid:12)(cid:12)(cid:1) , for α = ± w ( τ ) = 2 π [ ν ( τ )] d [ ν ( τ )] d + [ ν (2 π − τ )] d , ≤ τ ≤ πν ( τ ) = (cid:18) d − (cid:19) (cid:18) π − τπ (cid:19) + 1 d (cid:18) τ − ππ (cid:19) + 12 (46)The derivatives of w ( τ ) vanish up to order d − at theendpoints, and therefore d − derivatives of ξ α ( τ ) also vanishat τ = 0 , corresponding to ξ α (0) = α . Now, since J p,a ( a = u, v ) is expanded in terms of Chebyshev polynomials, whichsatisfy a discrete orthogonality property on the Chebyshevgrid points, we can accurately precompute the action of the K qpba , MT s,qpba and MT h,qpba operators on each Chebyshevpolynomial individually: K qpba [ T mn ] ( u (cid:48) , v (cid:48) ) = N βv − (cid:88) k =0 N βu − (cid:88) l =0 A qpba ( u (cid:48) , v (cid:48) , ξ u (cid:48) ( s l ) , ξ v (cid:48) ( t k )) ∂u∂s ( s l ) ∂v∂t ( t k ) w l w k T mn ( ξ u (cid:48) ( s l ) , ξ v (cid:48) ( t k )) (47)where ∂u∂s → and ∂v∂t → as ξ u (cid:48) ( s ) → u (cid:48) and ξ v (cid:48) ( t ) → v (cid:48) respectively, canceling the singularity in A . Note that theexpressions for MT s,qpba and MT h,qpba are the same but with A replaced by B and C respectively. It is important that N u,vβ is chosen sufficiently large to accurately compute each of theprecomputation integrals in (47) above. A numerical analysisof the resulting forward map accuracy vs N u,vβ is done inSection IV. Finally, on the basis of these precomputations, theaction of each of these operators on any J p,a or M p,a canbe readily computed using the Chebyshev expansion of thedensity, eg. K qpba [ J p,a ] ( u (cid:48) , v (cid:48) ) = N pv − (cid:88) m =0 N pu − (cid:88) n =0 γ p,am,n K qpba [ T mn ] ( u (cid:48) , v (cid:48) ) (48)where γ p,am,n are the Chebyshev expansion coefficients definedin (19). An analogous relation also holds true for the MT s and MT h operators. This precomputation approach is alsoused in order to accurately compute the K qpba , M T s,qpba and
M T h,qpba blocks corresponding to target points which are ondifferent patches but which are still in close proximity to thesource patch, making the integration near-singular. The onlydifference in this scenario arises in the selection of α in thechange of variable expression (46). Instead of simply choosingthe ( u (cid:48) , v (cid:48) ) corresponding to the target point, since it is on adifferent patch, we search for: ( u ∗ , v ∗ ) = arg min ( u,v ) ∈ [ − , | r q ( u (cid:48) , v (cid:48) ) − r p ( u, v ) | (49)for the change-of-variables as the point on the source patchnearest to the target patch, which can be readily found by anappropriate minimization algorithm. We adopted the goldensection search algorithm in our specific implementation [17]. IV. N UMERICAL R ESULTS
We first present the convergence of the forward map forboth the MFIE and N-M¨uller formulations with respect to theorder of expansion used ( N ) for varying levels of singularintegration refinement ( N β ). Following this, several numericalexamples involving scattering from PEC and dielectric spheresand cubes are presented and compared against a commercialRWG-based MoM solver to demonstrate the high accuracythat can be achieved using the proposed CBIE method. Fi-nally, we present scattering and near-field density results fromscattering by highly intricate 3D NURBS objects parametrizedwith commercial CAD software [18], which shows that theapproach can be readily applied to simulate objects arising inrealistic applications. A. Forward Map Convergence (a)(b)
Fig. 3: (a) Forward mapping error with respect to N forvarious choices of N β on a PEC sphere ( D = 2 λ ) using theMFIE formulation. (b) Forward mapping error on a dielectricsphere ( D = 2 λ e , (cid:15) e = 1 . , (cid:15) d = 2 . ) using the N-M¨ullerformulation.Fig. 3 plots the forward mapping error (i.e., applying theintegral operators to a reference solution and checking theaccuracy of the result) on a λ ( λ e ) diameter sphere geometryfor both the PEC and dielectric cases versus N for various different choices of N β . In the dielectric case, the exterior (cid:15) e = 1 . and the interior (cid:15) d = 2 . . The Mie series solution dueto an incident plane wave is used as the reference solution [19].As can be seen, depending on the desired accuracy, it isimportant to choose N β judiciously such that it does notlimit the overall solution accuracy. Increasing N β does notincrease the number of unknowns (controlled by N ); however,it can significantly increase the amount of time required toprecompute the singular and near-singular interactions. B. PEC scattering: MFIE Formulation
In this subsection, we test the proposed approach for theMFIE formulation by computing scattered fields from threePEC objects: two spheres of diameters . λ and λ and acube with side length . λ . All three objects are parameterizedby using 6 patches, and each patch is discretized with thesame number of unknowns N = N u = N v . Thus the totalnumber of unknowns per problem is Q = 2 × × N .The spheres are illuminated by the same plane wave source, E inc = exp ( − ikz ) ˆx . Since a closed-form solution does notexist for scattering from a cube, we use an electric dipoleexcitation, H inc ( r ) = −∇ × { G ( r , r (cid:48) ) p } , placed at position r (cid:48) = (0 . λ, . λ, . λ ) inside the cube with polarization p = (1 , , . This allows us to determine convergenceof the numerical solution since the scattered electric fieldmust cancel the incident field outside the cube, and thus: H scat ( r ) = ∇ × { G ( r , r (cid:48) ) p } (cid:0) r ∈ R \ D (cid:1) . The results for thesphere cases are compared against the analytical Mie seriessolutions.Fig. 4(a) shows the computed surface current density mag-nitude on the λ sphere for N = 26 and Fig. 4(b) shows theerror difference in surface density between the computed andanalytical solution on the λ sphere for N = 26 . As can beseen, the numerical solution differs from the exact solution byless than . × − at every point on the sphere. Fig. 4(c) plotsthe computed surface current distribution on the cube resultingfrom the internal dipole source. Fig. 5 plots the computedRCS overlaid on top of the analytical far-field solution of the λ sphere for varying altitudes at φ = 0 ◦ and φ = 90 ◦ . Thecomputed RCS points are indistinguishable from the analyticalsolution for both curves.Fig. 6 plots the error of the CBIE method vs number of un-knowns ( Q ) used to discretize each scatterer. As a comparison,the convergence of a commercial MoM RWG-based solver forthe λ sphere case is also plotted. For reference, st and th order slopes are drawn in dashed lines. As can be seen, theMoM solver only approaches first order convergence, requiresa much finer discretization than the proposed CBIE method,and even for a very high resolution mesh barely exceedstwo digits of accuracy. In contrast, CBIE converges spectrallyfast for all 3 examples, which makes it a significantly moreaccurate and efficient approach. C. Dielectric scattering: N-M¨uller formulation
The scattered fields from two dielectric objects are com-puted to evaluate the performance of the CBIE method for theN-M¨uller formulation: a dielectric sphere of 2 λ e diameter with (a) (b)(c) Fig. 4: (a) Surface current distribution on λ diameter sphere.(b) Error in surface current distribution on λ diameter sphere.The worst error is . × − . (c) Surface current distributionon . λ edge length cube.Fig. 5: Numerical vs exact RCSs corresponding to scatteringof plane wave from D = 4 λ PEC sphere at φ = 0 ◦ and φ = 90 ◦ .permittivity (cid:15) d = 2 (cid:15) e and a dielectric cube of 2 λ e side lengthwith permittivity (cid:15) d = 2 (cid:15) e , where the λ e = πk e is the wave-length corresponding the background exterior medium whichis set to free-space for all problems considered here ( (cid:15) e = (cid:15) ).The permeability for both objects is also set to the vaccumpermeability: µ d = µ e = µ . The surfaces of the objects arediscretized in the same manner as for the MFIE formulation,which results in Q = 2 × × × N unknowns. They are bothilluminated by a plane wave excitation E inc = exp ( − ikz ) ˆx .The results are compared against the Mie series analytical Fig. 6: Convergence of far-field error for the three scattererexamples vs number of unknowns. Performance of commercialMoM RWG-based solver is also plotted for D = 4 λ spherecase for comparison. nd and th order asymptotes are drawnfor reference. (a) (b)(c) (d) Fig. 7: (a) Surface J distribution on λ e diameter dielectricsphere with (cid:15) d = 2 (cid:15) e . (b) Error of surface J distribution. Maxerror: . × − . (c) Surface M distribution. (d) Error ofsurface M distribution. Max error: . × − .solution for the dielectric sphere [19] and against a highlyrefined numerical solution for the dielectric cube since ananalytical solution does not exist.Fig. 7(a) shows the electric (J) current density distributionon the surface of the λ e sphere for N = 24 and Fig. 7(b)shows the error difference of the computed electric currentdensity distribution with the Mie Series solution. Fig. 7(c) andFig. 7(d) show the magnetic (M) current density distribution and the corresponding error distribution respectively. It can beseen that the maximum relative error is smaller than . × − among all the discretization points on the sphere. Fig. 8 plotsthe computed RCS overlaid on top of the analytical far-fieldsolution of the λ sphere for varying altitudes at φ = 0 ◦ and φ = 90 ◦ .Fig. 9 plots the error of the CBIE method vs number ofunknowns ( Q ) used to discretize each scatterer. As expected,the convergence for the cube is considerably worse than that ofthe sphere due to the edge and corner singularities which causethe surface densities to not be smooth. The convergence ratecan be recovered, however, by using the same edge refinementapproach proposed in [9] which clusters unknowns near theedges to better resolve the singularities. This improvementcan be seen in the edge refined curve plotted in Fig. 9. Asa comparison, the convergence of a commercial MoM RWG-based solver for the both objects is also plotted. For reference, st and th order slopes are drawn in dashed lines. As withthe PEC case, the MoM solver only approaches first orderconvergence and requires a much finer discretization than theproposed CBIE method due to the linear basis functions andflat triangular discretization used to represent the geometry.Fig. 8: Numerical vs exact RCSs corresponding to scatteringof plane wave from D = 2 λ e dielectric sphere at φ = 0 ◦ and φ = 90 ◦ . D. Scattering from complex NURBS CAD models
In order to demonstrate that the proposed approach can bereadily used to solve scattering from complex CAD generatedmodels with arbitrary curvature, we solve for the scatteredfields from two different NURBS models freely available fordownload online [20]. As in the previous examples, the inci-dent excitation is an x -polarized plane wave propagating in the + z direction. In the first example, we consider scattering offof a 16 wavelength tall humanoid bunny character. Fig. 10(a)shows the induced surface current density and Fig. 10(b) plotsthe RCS vs θ at φ = 90 ◦ angle for two different discretizations( N = 10 and N = 12 Chebyshev points per side per patch or100 and 144 points per patch total respectively). The modelis comprised of 402 curvilinear quadrilateral patches total
Fig. 9: Convergence of far-field error for the two dielectricscatterer examples vs number of unknowns. Convergence forthe dielectric cube using edge refinement is also plotted.Performance of commercial MoM RWG-based solver is shownfor comparison. nd and th order asymptotes are drawn forreference. (a) (b) Fig. 10: (a) Surface electric current density induced on λ tall PEC CAD humanoid bunny model by incident plane wave.The model consists of 402 curvilinear quadrilateral NURBS-parametrized patches. (b) RCS at φ = 90 ◦ corresponding toplane wave scattering for N = 10 and N = 12 Chebyshevpoints per patch discretizations.and was directly imported from a standard CAD softwarewithout any special post-processing required [18]. Despite thelarge size of the model, significant variation in curvature, andregions with sharp corners (eg., the ears), the match in the RCSfor the two relatively coarse discretizations is excellent andthey are almost indistinguishable from one another, varyingless than × − from each other.For the second CAD model example, we computed scat-tering from a glider with a length of 7.7 wavelengths and awingspan of 5.6 wavelengths from the end of one wing to theother. Fig. 11(a) shows the induced surface current densityand Fig. 11(b) plots the RCS vs θ at φ = 90 ◦ angle for two different discretizations ( N = 10 and N = 12 Chebyshevpoints per side as before). The glider is comprised of 79curvilinear quadrilateral patches total. As before, the RCScurves resulting from the two different discretizations matchvery well and vary less than . × − from each other. (a)(b) Fig. 11: (a) Surface electric current density induced on 79patch PEC glider CAD model by incident plane wave. Theglider spans 8 wavelengths from wing to wing. (b) RCS at φ = 90 ◦ corresponding to plane wave scattering for N = 10 and N = 12 Chebyshev points per patch discretizations.V. C
ONCLUSION
This paper presents a high-order accurate Chebyshev-basedBoundary Integral Equation (CBIE) approach for solvingMaxwell’s equations. The CBIE method is applied towards thediscretization of the MFIE and the N-M¨uller formulation. Theperformance is evaluated by solving scattering from sphere andcube PEC/dielectric objects and comparing against analyticalsolutions as well as a commercial MoM-based solver. We havealso demonstrated a couple examples of scattering from com-plex 3D CAD models which contain many intricate featuresand variations in curvature. The proposed method achievesspectral convergence on sufficiently smooth surfaces withrespect to the number of unknowns, significantly reducing thenumber of unknowns required for a desired accuracy over low-order MoM approaches. Furthermore, the CBIE approach alsoconverges well for geometries with edges and corners when an edge-refinement change of variables is utilized as demon-strated by the dielectric cube example. Current and futurework involves applying the CBIE method in conjunction withthe Windowed Green Function (WGF) [21] method towardsthe simulation and design of 3D waveguiding structures withunbounded boundaries for modeling nanophotonic devices [5],treating multi-material and composite objects [22]–[24], andincorporating acceleration techniques such as the Fast MultipleMethod [25] or FFT-based methods [26], [27].R EFERENCES[1] J. J. Bowman, T. B. Senior, and P. L. Uslenghi, “Electromagnetic andacoustic scattering by simple shapes (revised edition),” in hpc , 1987.[2] C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch,“Adjoint shape optimization applied to electromagnetic design,”
Opticsexpress , vol. 21, no. 18, pp. 21 693–21 701, 2013.[3] A. Taflove and S. C. Hagness,
Computational electrodynamics: the finite-difference time-domain method . Artech house, 2005.[4] O. C. Zienkiewicz, R. L. Taylor, P. Nithiarasu, and J. Zhu,
The finiteelement method . McGraw-hill London, 1977, vol. 3.[5] C. Sideris, E. Garza, and O. P. Bruno, “Ultrafast simulation andoptimization of nanophotonic devices with integral equation methods,”
ACS Photonics , vol. 6, no. 12, pp. 3233–3240, 2019.[6] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering bysurfaces of arbitrary shape,”
IEEE Transactions on antennas and prop-agation , vol. 30, no. 3, pp. 409–418, 1982.[7] F. P. Andriulli, “Loop-star and loop-tree decompositions: Analysis andefficient algorithms,”
IEEE Transactions on Antennas and Propagation ,vol. 60, no. 5, pp. 2347–2356, 2012.[8] O. Bruno, T. Elling, R. Paffenroth, and C. Turc, “Electromagnetic inte-gral equations requiring small numbers of krylov-subspace iterations,”
Journal of Computational Physics , vol. 228, no. 17, pp. 6169–6183,2009.[9] O. P. Bruno and E. Garza, “A chebyshev-based rectangular-polar integralsolver for scattering by general geometries described by non-overlappingpatches,” arXiv preprint arXiv:1807.01813 , 2018.[10] A. Maue, “On the formulation of a general scattering problem by meansof an integral equation,”
Z. Phys , vol. 126, no. 7, pp. 601–618, 1949.[11] C. M¨uller,
Foundations of the mathematical theory of electromagneticwaves . Springer Science & Business Media, 2013, vol. 155.[12] P. Yl¨a-Oijala, M. Taskinen, and S. J¨arvenp¨a¨a, “Analysis of surfaceintegral equations in electromagnetic scattering and radiation problems,”
Engineering Analysis with Boundary Elements , vol. 32, no. 3, pp. 196–209, 2008.[13] J. Volakis,
Integral equation methods for electromagnetics . TheInstitution of Engineering and Technology, 2012.[14] P. Yla-Oijala and M. Taskinen, “Well-conditioned muller formulationfor electromagnetic scattering by dielectric objects,”
IEEE Transactionson Antennas and Propagation , vol. 53, no. 10, pp. 3316–3323, 2005.[15] J. C. Mason and D. C. Handscomb,
Chebyshev polynomials . CRCpress, 2002.[16] D. L. Colton, R. Kress, and R. Kress,
Inverse acoustic and electromag-netic scattering theory . Springer, 1998, vol. 93.[17] J. Kiefer, “Sequential minimax search for a maximum,”
Proceedings ofthe American mathematical society
IEEETransactions on Antennas and Propagation , vol. 65, no. 9, pp. 4684–4692, 2017.[22] P. Yla-Oijala, M. Taskinen, and J. Sarvas, “Surface integral equationmethod for general composite metallic and dielectric structures withjunctions,”
Progress In Electromagnetics Research , vol. 52, pp. 81–108,2005.[23] C. P´erez-Arancibia, C. Turc, L. Faria, and C. Sideris, “Planewave densityinterpolation methods for the efie on simple and composite surfaces,” arXiv preprint arXiv:1910.02046 , 2019.[24] C. P´erez-Arancibia, C. Turc, and L. Faria, “Planewave density inter-polation methods for 3d helmholtz boundary integral equations,”
SIAMJournal on Scientific Computing , vol. 41, no. 4, pp. A2088–A2116, 2019. [25] N. Engheta, W. D. Murphy, V. Rokhlin, and M. S. Vassiliou, “The fastmultipole method (fmm) for electromagnetic scattering problems,”
IEEETransactions on Antennas and Propagation , vol. 40, no. 6, pp. 634–641,1992.[26] O. P. Bruno and L. A. Kunyansky, “A fast, high-order algorithm for thesolution of surface scattering problems: basic implementation, tests, andapplications,”
Journal of Computational Physics , vol. 169, no. 1, pp.80–110, 2001.[27] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, “Aim: Adaptiveintegral method for solving large-scale electromagnetic scattering andradiation problems,”