A Concise and Efficient MLFMA Scheme for Electromagnetic Surface Integral Equations From Dielectric Objects

In the surface field integral equations (SIEs) from three-dimensional homogeneous dielectric objects, the equivalence electric and magnetic currents on the discontinuous interfaces between different homogeneous media are two independent sources. When the SIEs are iteratively solved, as indicated by the published literatures, the multilevel fast multipole algorithm (MLFMA) needs to be implemented respectively for the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) with these two currents, and hence a few MLFMA implementations are necessary for each homogeneous medium. In this paper, a concise and efficient scheme is proposed for reducing the implementations of the SIEs-MLFMA. First, the summation of the EFIE and MFIE-MLFMA implementations with the two currents can be simplified into one EFIE or MFIE-MLFMA implementation with integrating electric and magnetic currents at the lowest-level aggregation stage. Secondly, the common matrix frameworks consisting of the EFIE and MFIE-MLFMA implementations used in the SIEs can be simplified into one MLFMA implementation with two lowest-level diaggregations. The enhanced MLFMA scheme has excellent applicability to eight kinds of SIEs in wide use. In this way, the resulted SIE-MLFMA scheme performs only once for each homogeneous medium when transforming outgoing wave expansions of the two currents into incoming wave expansions for bottom cubes. Compared with the widely-used implementation approach of the MLFMA one by one for the two currents, it reduces the computation complexity by one half without sacrificing accuracy and adding a little memory requirement. Some examples are presented to validate the enhanced MLFMA scheme.


I. INTRODUCTION
Surface integral equations (SIEs) have been widely used for formulating electromagnetic (EM) analysis of objects made of metallic and piecewise homogeneous dielectric media [1], [3] since the equivalence electric and magnetic current sources J and M distribute only on the discontinuous interfaces between different media, which leads to less unknowns than the volume integral equations. The tangential EFIE and MFIE equations and their rotational versions n × EFIE and n×MFIE are basis SIEs for each homogenous The associate editor coordinating the review of this manuscript and approving it for publication was Su Yan . medium. Their integral operators L and K, and rotL and rotK constitute two matrix frameworks, i.e., the 2 × 2 matrix of L and K, and the matrix of rotL and rotK. Weighting these matrics with different scales yields a variety of coupled SIE formulations [3]- [5], [7], [8], [14], [21] such as PMCHWT, Müller, CTF, CNF, and JMCFIE, which are free of the internal resonance problem. The PMCHWT formulation [3], [21] yields an accurate solution with a slow convergence rate. The N-Müller formulation [8], [9] can achieve high accuracy in the low-frequency case but results in large error in the high-frequency case, in particular, for objects with high permittivity and sharp edges. It converges much faster than the PMCHWT. By normalizing field quantities and electric and VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ magnetic currents [3], the JMCFIE and its variant versions with scaled matrices such as TSJMCFIE, RSJMCFIE, and SSJMCFIE formulations developed in [3] work well on broad ranges of permittivities and frequencies, and have fast convergence rates with sacrificing a little accuracy compared with the PMCHWT and CTF formulations [4]. But their solutions are more accurate than the CNF and N-Müller formulations which are contaminated with the inaccuracies due to the low-order discretizations of the identity operator [19], [20]. Hence, a tradeoff between the efficiency and accuracy needs to be evaluated when the SIEs are used for formulating EM problems from dielectric objects. There are two unknowns on each RWG element [10] corresponding to J and M. As a result, the MoM matrix equations of the discretized SIEs are memory-consuming and the iteration solutions are time-consuming, which would be the bottleneck for large-scale EM problems. Various fast methods such as the fast Fourier transform (FFT) [11], the adaptive integral method (AIM) [12] and the multilevel fast multipole algorithm (MLFMA) [6], [8], [13]- [17], [22]- [29] have been developed. The MLFMA is one of the great breakthroughs of computational electromagnetics (CEM), and it reduces the computation time and memory requirement to O(NlogN ) for a matrix-vector multiplication (MVM) [13], where N is the number of degrees of freedom.
When the MLFMAs are applied to the coupled tangential EFIE and MFIE formulations such as PMCHWT and CTF, which have the first matrix framework, four MLFMA implementations of L and K with independent J and M are needed for each homogeneous medium, which are fourfold of that in the case of PEC objects with same size. Hence, the MLFMA in the case of dielectric objects is much less efficient than that in the PEC case. In [4], [5], [13], and [14], the first column in the first matrix framework, i.e., a pair of MLFMA implementations of L and K operators with J (or the second column, i.e., a pair of L and −K with M) are integrated into one MLFMA implementation for each homogeneous medium but with two lowest-level disaggregations of incoming wave expansions from the bottom cube to the testing RWG functions. In the case of rotational EFIE and MFIE formulations such as N-Müller and CNF, which have the second matrix framework, the corresponding first and second columns can be integrated into two MLFMA implementations, too.
In the combined versions of tangential and rotational EFIEs and MFIEs such as JMCFIE, TSJMCFIE, SSJMCFIE, and RSJMCFIE [3]- [5], which have the two matrix frameworks, two MLFMA implementations with four lowest-level disaggregations are also needed in each homogeneous medium. However, such two MLFMA implementations are still expensive especially when analyzing electrically-large problems. An efficient version would still be desirable, which motivates our work in this paper.
Based on the vector operation identities, we simplify the first matrix framework consisting of L and K operators with J and M into one column consisting of L(J +k × M) and K(J +k × M) from view point of MLFMA, and the second matrix framework of rotL and rotK with J and M into one column of rotL(J +k × M) and rotK(J +k × M). Obviously, the two column frameworks have the same aggregations, translations, disaggregations from the top level to the lowest level. At the lowest level, their differences exist in the disaggregations of wave expansions from the bottom cube to the RWG functions, the n× RWG functions, thek× RWG functions, and thek × n× RWG functions. Hence, in the column frameworks, the MLFMA is implemented only once at the aggregation, translation, and disaggregation stages except for the lowest-level disaggregation.
Due to containing either one or two common column frameworks, the MLFMA schemes of the JMCFIE, TSJM-CFIE, SSJMCFIE, RSJMCFIE, PMCHWT, N-Müller, CTF, and CNF can be expressed as concise and efficient MLFMA schemes, respectively. This paper is organized as follows. Section II describes the normalized CTF and CNF and their matrix frameworks, which constitute various SIEs. In Section III, by integrating electric and magnetic currents at the lowest-level aggregation stage, the two matrix frameworks consisting of the MLFMA implementations of the tangential and rotational EFIEs and MFIEs are converted into an efficient and concise MLFMA scheme with a few lowest-level disaggregations for each homogeneous medium. In Section IV, the enhanced MLFMA scheme is applicable to the PMCHWT and N-Müller formulations, respectively. Section V outlines the various JMCFIE formulations accelerated by the enhanced MLFMA scheme. In Section VI, some numerical examples are presented to demonstrate the capability of the presented scheme. Finally, several conclusions are given in Section VII.

II. FORMULATIONS
Consider a dielectric object residing in a homogeneous background medium. The outer and inner media of the object are characterized with EM parameters i and µ i for i = 1, 2, respectively. The EM problem can be formulated by the CTF and CNF [3] in terms of the normalized equivalence surface electric and magnetic currents where n i is the unit normal vector pointing into the outer and inner (i = 1, 2) media on the equivalence surface S of the dielectric object. The CTF and CNF formulations in the outer and inner media can be expressed, respectively, as where rotL i and rotK i denote n i ×L i and n i ×K i , respectively. Either one or two frameworks are shared by the SIEs such as the PMCHWT, N-Müller, CTF, CNF, and various JMCFIE formulations. The electric and magnetic currents are expanded with the RWG functions [10], respectively, Testing the PMCHWT formulation with the RWG functions yields the impedance matrix equations where the near-field matrix Z PM near represents the near-field interactions of the electric and magnetic currents in the neighboring bottom cubes in the MLFMA, respectively. It sums up the near-field interactions from the inner and outer media simultaneously, and its elements are calculated through the approach in [2]. The far-field matrices represent the far-field interactions of the current sources in the well-separated cubes [13]. For i = 1, 2, L i and K i denote the far-field submatrices of L i and K i in the matrix framework in (9), respectively. They are not numerically available and implicit in the MLFMA.

the excitation vectors of E inc
1 and H inc 1 tested with the RWG functions. The discretized N-Müller can be written as where rotL i and rotK i denote the far-field submatrices of rotL i and rotK i for i = 1, 2 in the matrix framework in (10), and rotb H and rotb E are similarly defined. During the process of the iteration solution to either (12) or (13), the near-field matrix-vector multiplications (MVMs) are directly implemented, and the far-field MVMs are accelerated by the MLFMA. Since the Green function G i relies on the wavenumber k i for i = 1, 2, the aggregation, translation, and disaggregation after factorizing and diagonalizing G i in the MLFMA are intensively dependent on k i [13]. Consequently, more than one MLFMA are needed to accelerate the far-field MVMs in (12) and (13) for the outer and inner homogeneous media, respectively.

III. MLFMA WITH ELECTRIC AND MAGNETIC CURRENT SOURCES
For the sake of clarity, the subscript ''1'' or ''2'' in (12) and (13) is omitted in this section. We pay our particular attention to the MVMs in the far-field matrix frameworks in (12) and (13) as following and The vectors used in these MVMs in (14) and (15) are still denoted by the notations in (12) and (13) for simplicity if such treatment does not lead to confusion. The uth element of Lx in (14) can be expressed as which d 2k denotes the integral over the unit sphere, and L the lowest level of the octo-tree. The incoming wave expansion from the center whichĪ is the unit dyadic andk = k/|k|; S u is the RWG element. When 2 ≤ l ≤ L, and where {N m l } denotes the set of the cubes neighboring to C m l . The incoming wave expansions for C m l in (18) consist of two parts: ones are disaggregated from the incoming wave expansions W l−1 for the parent cube C m l−1 with the anterpolation technique, and the others are translated from the outgoing wave expansions V l,source n l for the non-neighboring cube C n l whose parent cube C n l−1 is neighboring to C m l−1 . The diagonal translation operator [13], [22] is denoted by α m l n l .
The outgoing wave expansion for C n l is where x v is the element of x in (14) and the surface integral on the right hand side represents the outgoing wave expansions for each RWG function f v in C n L are aggregated to the center of C n L . The outgoing wave expansions for the cubes at low level are gradually aggregated to those for the cube at high level by the interpolation technique. (20) is replaced with y v of y, the uth element of Ly in (14) can be expressed as Sincek the uth elements of Kx and Ky in (14) can be written, respectively, as and A. MLFMA RESPECTIVELY WITH ELECTRIC AND MAGNETIC CURRENTS Concluded from (16), (21), (23) and (24), the MLFMA schemes for both the uth and (N + u)th elements of MVMs in (14) can be written as where the incoming wave expansions for C m L , i.e., W L (k, J) or W L (k, M), are simultaneously disaggregated to the uth RWG function with different testing choices. Similarly, the MLFMA schemes for both the uth and (N + u)th elements of MVMs in (15) can be written as As suggested in [5], [13], [14], W L (k, J) in (25) and (26) is preferably multiplied with the first columns of the matrices consisting of different V L,field m L , and W L (k, M) multiplied with the second columns. In such way, integrating the MLFMA implementations into columns is feasible.
Hence, two aggregations, two translations, and two disaggregations from level-2 to level-L are necessary due to independent sources J and M. Moreover, four lowest-level disaggrations of incoming wave expansions from the bottom cube C m L to f u are essential due to different testing procedures. Roughly speaking, two MLFMA implementations are required for accelerating far-field MVMs in (14) or (15) for each homogeneous medium, which can be further reduced in what follows.

B. MLFMA INTEGRATING ELECTRIC AND MAGNETIC CURRENTS
Deduced from (22) and (18)- (20), the MLFMA scheme (24) for the uth element of Ky can be transformed to which means that ''k×'' is moved equivalently from the lowest-level disaggregation of incoming wave expansion into the lowest-level aggregation of outgoing wave expansion. It is noted that thek direction samples over the unit sphere at the level-l aggregation and disaggregation are the same for l = 2, · · · , L. Because of common V L,field m L in (16) and (27), the uth element of [Lx − Ky] in (14) can be expressed as where Having similar form to (16), (28) may represent the far-field interactions of integrated currents J +k × M on f u for L operator, intuitively.
Sincek × f u ×k = (Ī −kk) · f u , the MLFMA scheme for the uth element of Ly can be transformed to Obviously,k × f u remains in the lowest-level disaggregation of incoming wave expansions, and ''k×'' is moved into the lowest-level aggregation of outgoing wave expansions. Because of common V L,field m L , summing (23) and (29) (14), as which represents the far-field interactions of J +k × M on f u for K operator, intuitively. Because of common W L , (28) and (30) can be further simplified into an efficient and concise MLFMA scheme for both the uth and (N + u)th elements of MVMs in (14), In (31), the MLFMA performs once from the aggregation of outgoing wave expansions of integrated currents J +k × M to the disaggregation of incoming wave expansions for C m L , and the two lowest-level disaggregations from C m L to f u are required due to different testing choices. In similar way, an efficient and concise MLFMA scheme for both the uth and (N + u)th elements of MVMs in (15) can be written as Obviously, the MLFMA schemes (31) and (32) have the same aggregation, translation, and disaggregation from level-2 to level-L, and their differences exist in the lowest-level disaggregation of incoming wave expansions from the bottom cube C m L to f u . This common property facilitates these MLFMA schemes to be integrated into a very concise form. The enhanced version performs with less computation complexity for accelerating the far-field MVMs from the discretized PMCHWT, N-Müller, and different JMCFIE versions [3] containing both (14) and (15).
Fundamentally, both (31) and (32) are successful cases of inverse application of distributive law of vector operations to aggregation, translation and disaggregation. However, both (25) and (26) are incomplete cases.

IV. PMCHWT-MLFMA AND N-MÜLLER-MLFMA
Because of the far-field matrix framework (14) contained in (12), utilizing (31), the outer MLFMA scheme for both the uth and (N + u)th elements of the first far-field MVMs in (12) can be expressed as The inner version for both the uth and (N + u)th elements of the second far-field MVMs in (12) is given by Summing (33) and (34) yields both the uth and (N + u)th elements of the far-field MVMs in (12). In such way, two MLFMAs are sufficient to accelerate the far-field interactions.
In the popular approach depicted in [4], [5], and [14], the outer MLFMA schemes with independent J 1 and M 1 for both the uth and (N + u)th elements of the first far-field MVMs in (12) are given by (25) with adding subscript ''1'', and the inner versions for both the uth and (N + u)th elements of the second far-field MVMs in (12) are given by Obviously, (33) and (34) reduce the computation complexity of the far-field MVMs by almost one half without additional memory requirement compared with (35) and (36).
Since the far-field matrix framework (15) is contained in the discretized N-Müller formulation (13), based on (32), the outer MLFMA scheme for both the uth and (N + u)th elements of the first far-field MVMs in (13) can be expressed as The inner version for both the uth and (N + u)th elements of the second far-field MVMs in (13) is given by Combining (7) and (8) yields the JMCFIE formulation with normalized field quantities. Its discretized matrix equations can be obtained. In particular, the far-field interactions read Summing (31) and (32) yields the outer and inner MLFMA schemes for both the uth and (N + u)th elements of the far-field MVMs in (39) as The concise TSJMCFIE-MLFMA scheme is given for the outer and inner media (i = 1, 2) as When s µ i and s i exchange places in (41), and are replaced with their inverses s −1 i and s −1 µ i , respectively, the enhanced SSJMCFIE-MLFMA and RSJMCFIE-MLFMA schemes can be gotten, respectively.
The discretized SIEs mentioned above have either one or two common far-field matrix frameworks, i.e., (14) or/and (15). Then, the highly-efficient MLFMA schemes (31)-(32) are available for building these SIEs-MLFMA schemes. In these schemes, only one MLFMA implementation with two lowest-level disaggregations is needed for each homogeneous medium, which shows that Sections III-V provide versatile electromagnetic SIE-MLFMA solvers with high efficiency.

VI. NUMERICAL EXAMPLES
Limited to length of content, as typical formulations, the TSJMCFIE-MLFMA scheme and the PMCHWT-MLFMA scheme integrating the electric and magnetic currents are tested on some dielectric objects through comparison with the Mie series [30] or those versions respectively with the electric and magnetic currents. Triangular patches with dimension λ/10 are used to mesh the surface of the  object, and λ denotes the wavelength in free space. All computations are carried out on a HP Z440 workstation with a 6-cores Intel Xeon processor and storage of 128 GB. In all cases, the incident wave is a plane wave with x-polarization and z-propagation. The scattering angle θ is measured from the z-axis in the x-o-z plane. The GMRES method [18] as an iteration solver is terminated when its relative residual error reaches the tolerance value of 10 −3 .
A line array consisting of 28 RWG elements is shown in Fig.1. The RWG elements are of the same size and are in the x-o-y plane in Cartesian coordinate system. The source RWG element is colored in light gray, and its center is located at the origin point. The centers of the field RWG elements colored in dark gray are located at x n = 0.35(n + 1)λ for n = 1, · · · , 27 along the x-axis. When the source RWG element is rotated along the x-axis, the angle α between the source and field RWG elements varies.
The size of bottom cube used here is set 0.35λ, and then the bottom cubes containing the field RWG elements are not neighboring to the bottom cube containing the source RWG element. Hence the interaction from the source RWG function on the 27 field RWG functions can be calculated by the the different MLFMA schemes in Section III. These calculations are under the conditions of quadrature rule, expansion, interpolation and anterpolation techniques.
The interactions between the source and field RWG functions for L operator are calculated by the EFIE-MoM, the MLFMA versions of (29) and (21), and their results and errors are compared respectively for α = 0 • , 45 • and 90 •  Since the cube containing the field RWG element located at x 1 = 0.7λ is well-separated [13] to the cube containing the source RWG element located at the origin point, without undergoing high-level aggregation and disaggregation stages, the outgoing wave expansions for the latter cube are directly translated to the former cube. Consequently, excellent agreement of results at x 1 from the MLFMAs schemes and the MoM can be observed from Figs. 2-5. However, tiny deviations of results exist at x n > 0.7λ. The reason is that the interpolation and anterpolation in the aggregation and disaggregation stages decrease the accuracy of (29) and (21) calculating L operator and (27) and (24)   operator, respectively. The other errors sources include the limited expansions of addition theorem of Green function, the Lagrange polynomial interpolation technique employed in the translation stage, and the quadrature integral along the θ and φ directions in d 2k .
A sphere of radius 5λ is firstly considered, where 226,662 unknowns are involved. Fig. 6 shows that the RCS results from the PMCHWT-MLFMA integrating J and M (i.e., (33) and (34)), and the PMCHWT-MLFMA respectively with J and M (i.e., (35) and (36)) are in well agreement with those from the Mie series [30], respectively. Fig. 7 compares the convergence performances of these two PMCHWT-MLFMA versions and shows that their iteration processes are of the same convergence rate. Seen from Fig. 7, these TSJMCFIE-MLFMA versions have the same convergence performances. Fig. 8 shows that the results from the TSJMCFIE-MLFMA integrating J and M (i.e., (41)) are identical to those from TSJMCFIE-MLFMA respectively with J and M. From point view of the iteration convergence history     A electrically-large airplane in real life is shown in Fig. 9, where 3,049,986 unknowns are involved. Fig. 10 shows that the RCS results from the MLFMA integrating J and M agree well with those from the MLFMA respectively with J and M for the TSJMCFIE. The PMCHWT-MLFMA fails to converge to the tolerance value of 10 −3 within 700 iterations. Seen from Fig. 11, both TSJMCFIE-MLFMA versions have the same fast-convergence performance. The CPU time of the iteration solution of the TSJMCFIE-MLFMA respectively with J and M is 294.78 minutes, and that of the TSJMCFIE-MLFMA integrating J and M is 154.57 minutes, which shows that the reduction of iteration time reaches 47.56%. Their memory requirements are 40.6 GB.

VII. CONCLUSION
In this paper, a concise and efficient MLFMA scheme is proposed for almost all kinds of SIEs involving homogeneous dielectric objects. Through vector operator identities, four individual MLFMA implementations of integral operators of tangential (or rotational) EFIE and MFIE with electric and magnetic currents are merged into single MLFMA implementation with integrated currents J +k × M for each homogeneous medium, which is capable of being seamlessly assembled with eight kinds of SIEs in wide use. The enhanced SIE-MLFMA scheme reduces the computational cost almost by one half without accuracy compromise and additional memory requirement, compared with the widely-used SIE-MLFMA version, and it can serve as a highly-efficient and concise approach to solving EM problems from large-scale dielectric objects. The enhanced SIE-MLFMA scheme can easily be extended to EM problems involving composite structures with multiple dielectric and metallic objects. Since 2001, she has been with the School of Applied Mathematics, Nanjing University of Finance and Economics. Her research interests include numerical methods in engineering problems, matrix analysis, and optimization method.
WEI HONG (Fellow, IEEE) received the B.S. degree from the University of Information Engineering, Zhengzhou, China, in 1982, and the M.S. and Ph.D. degrees from Southeast University, Nanjing, China, in 1985 and 1988, respectively, all in radio engineering.
Since 1988, he has been with the State Key Laboratory of Millimeter Waves and serves for the Director of the Lab, since 2003. In 1993In , 1995In , 1996In , 1997In , and 1998, he was a short-term Visiting Scholar with the University of California at Berkeley and at Santa Cruz, respectively. He is currently a Professor with the School of Information Science and Engineering, Southeast University. He has been involved in numerical methods for electromagnetic problems, millimeter wave theory and technology, antennas, RF technology for wireless communications, and so on. He has authored and coauthored over 300 technical publications and two books. Dr