Three-Dimensional Scattering From Uniaxial Objects With a Smooth Boundary Using a Multiple Infinitesimal Dipole Method

The formulations for three-dimensional (3D) scattering from uniaxial objects with a smooth boundary using a multiple infinitesimal dipole method (MIDM) are introduced. The proposed technique uses two sets of infinitesimal dipole triplets (IDTs), including three co-located orthogonally polarized electric infinitesimal dipoles, distributed inside and outside of a scatterer to construct simulated fields. The dyadic Green’s functions of uniaxial materials are deployed in the MIDM so as to obtain the simulated fields. The singularity issues in using the uniaxial dyadic Green’s functions, which cannot be solved analytically so far for a general uniaxial medium, can be easily eliminated by using the proposed MIDM. In comparison to the traditional single-layered distribution scheme of IDTs, the proposed multiple-layered distribution scheme can handle the scattering from uniaxial objects accurately and efficiently. Several numerical examples are presented to study bistatic radar cross section (RCS) responses under different scenarios. Excellent agreement is achieved by comparing numerical results with those obtained from commercial software packages, while the simulation performance including CPU time and required memory is drastically improved by using the MIDM when computing a general uniaxial material or a relatively larger object. The proposed technique has its merits on simplicity, conciseness and fast computation in comparison to existing numerical methods.


I. INTRODUCTION
The interaction between electromagnetic waves and anisotropic materials has received a great deal of attentions recently. Anisotropic materials have found a variety of applications in the design of antennas [1]- [7], integratedcircuit structures [8], reduction of RCS of scatterers [9], optical signal processing [10] and so on. One of the basic problems to investigate waves in the anisotropic material is The associate editor coordinating the review of this manuscript and approving it for publication was Su Yan .
to study the scattering characteristics of an anisotropic object. Several remarkable research contributions have been made and presented in [11]- [23]. The volumetric integral equation (VIE)-based methods were introduced in [11]- [16], [23] to compute scattering performances from an arbitrarily shaped object made of a linear, lossy, and anisotropic material. The VIE-based approaches can handle the most general cases of materials whereas they require to discretize the entire volume of an object, therefore becoming computationally challenging with large scatterers. The same issue will also rise in the finite-difference time domain (FDTD) [19], [22] and finite element-boundary integral (FE-BI) [20] methods. The surface integral equation (SIE) is a good candidate to overcome the computational burden of the methods based on volumetric discretization. A SIE-based MoM scheme combined with uniaxial dyadic Green's functions [24], [25] was proposed in [17], [18] for scattering evaluation from arbitrarily shaped objects filled with electrically uniaxial materials. The SIE-based solutions provide an accurate and more simplified approach compared to the VIE-based ones, yet the complexity of formulations is still there. The proposed analytical approach in [17] for eliminating singularity issues is only valid for an electrically uniaxial material, and it would fail when a general uniaxial material is encountered.
The generalized multipole technique (GMT) [26], [27] is a generic name of several similar numerical methods [28]- [31] developed independently by several research groups. In the GMT, the scattered fields are usually expanded in terms of a set of multipole sources. However, not only the multipoles can be used for fields expansion, but other equivalent sources are also possible. Therefore, other names for similar methodologies have been given like multiple multipole method (MMP) [28], discrete sources method (DSM) [29], method of auxiliary sources (MAS) [30], or multifilament current model (MFCM) [31]. A common basic concept of all these methods is that the scattered fields inside and outside of a scatterer are simulated by a set of equivalent sources respectively located outside and inside of the scatterer with a certain distance away from the physical boundary, rather than being formulated in terms of equivalent surface currents flowing on the physical surface. In this case, no integrals have to be computed numerically which reduces the computation time and simplifies the problem formulation. Also the solution features no singularity. So far, most of the scattering problems tackled with the GMT have considered isotropic objects, and there is few research work on the anisotropic scenario. In [32], the GMT was extended to anisotropic scatterers by introducing the plane wave representation of an anisotropic material into Bessel multipoles, but it resulted in integrals which cannot be evaluated analytically to represent the scattered fields. The DSM was extended to three-dimensional anisotropic scatterers in [33], but the entire body of the scatterer needed to be discretized. As a result, DSM also suffers from a high computational burden as in the case of the VIE-based method when larger objects are involved. Recently, the random auxiliary sources (RAS) method, a MMP-inspired numerical technique, was introduced in [34]- [36]. But all applications of the RAS currently only focused on isotropic materials. To the best of the authors knowledge, there is no research reported to date about combining the GMT-like methods with dyadic Green's functions of anisotropic materials to study the scattering from anisotropic objects.
The uniaxial material seems to be the most widely used type of anisotropic materials. This is because the uniaxial material can be either easily found in many natural crystals [8], [37], or artificially made by a stacked dielectric sheet structure consisting of alternative layers of two isotropic materials [4], [5], [38], [39], or obtained by homogenizing a mixture of several different materials via effective medium theory [40], [41]. This work only focuses on the uniaxial material, yet the proposed MIDM can also be applied to other kinds of anisotropic materials as long as the corresponding dyadic Green's functions are available. In the MIDM, the integral operations used in the VIE-or SIE-based methods are avoided. Instead, a simple point-matching testing procedure is used and the harmonics of a multipole source which is used in the MMP are also cast off by deploying the infinitesimal dipole source. Without using the integral operations and the harmonics of the sources, the computation of scattered fields speeds up and the problem formulation is also simplified. However, in comparison with the MMP which uses high order harmonics in equivalent sources, a larger number of infinitesimal dipole sources will be required in the MIDM to represent varying fields. An approach similar to the proposed MIDM is used in [42]. However, only isotropic materials were considered, and the strategy for the placement of sources may fail when a relatively larger object is involved, as will be shown in Section IV.
The scattering from uniaxial objects with a smooth boundary using a MIDM is proposed in this paper. The applications of structures with a smooth boundary can be found in many scenarios, such as the spherical dielectric resonator antennas, lens antennas, extended hemispherical lens antennas and so on. The paper is organized as follow. The formulation of the problem is presented in Section II, where the strategies for placements of matching points and sources that play a key role in the MIDM, will be discussed in detail and specified there. The singularity issues in using uniaxial dyadic Green's functions are discussed and solved in Section III. Several numerical examples are presented in Section IV under different scenarios. All numerical results are compared with simulated results generated from commercial software package, and excellent agreement is obtained. Finally, a conclusion to summarize the proposed technique is given in Section V.

II. FORMULATION OF THE MIDM
The geometry of the problem is illustrated in Fig. 1. Two regions are present, the external region 1 is free space and The rotated global coordinate system, represented by unit vectorsû,v andĉ as shown in Fig. 2, are used to express the uniaxial medium.ĉ is the unit vector parallel to the optical axis of the medium. 1 , µ 1 and 2 , µ 2 are the permittivity and permeability associated with the directions perpendicular and parallel to the optical axis (ĉ), respectively. The relationship between the unit vectors in the rotated global and global coordinates is written as  where θ c and ϕ c are defined in Fig. 2. The concept of the MIDM is illustrated in Fig. 3. We place a set of infinitesimal dipole triplets (IDTs) in regions 1 and 2. Each IDT contains three co-located orthogonal polarized infinitesimal dipoles, namely point sources. The formulation of the proposed MIDM is conducted through a two-step equivalence. Firstly, the scattered fields in the region 1 are generated by equivalent IDTs placed in the region 2, and those point sources are treated as source currents radiating in unbounded vacuum. Secondly, the internal fields in the region 2 are generated by equivalent IDTs placed in the region 1, and those point sources are radiating in unbounded space filled with an homogeneous uniaxial material identical to that constituting the scatterer.

A. FIELDS EXPRESSIONS IN REGIONS 1 AND 2
Region 1, considered as free space, contains the incident (E inc , H inc ) and the scattered fields (E s 1 , H s 1 ). If the incident fields are the plane wave, they could be written as: E inc =pe −jk(x sin θ inc cos ϕ inc +y sin θ inc sin ϕ inc +z cos θ inc ) (3a) wherep is the polarized direction of the incident electric field, andk is the unit vector of the wave vector. The scattered fields could be constructed in a simple manner and expressed in terms of the dyadic Green's functions. So the total fields (E 1 , H 1 ) could be expressed as: where J 1 i1 , J 1 i2 and J 1 i3 are the three orthogonal electric point sources in the ith IDT associated to the region 1. Noticing that three magnetic point sources can also be deployed in each IDT, and dyadic Green's functions (G em and G mm ) are then required with respect to the magnetic point source. Only electric point source based IDT is utilized throughout this paper. In (4), N 1 is the number of IDTs placed in the region 2. The two Green's functions, G 1 ee and G 1 me , in (4) are the dyadic Green's functions of isotropic materials, corresponding to the electric and magnetic fields radiated into the region 1 by an electric point source, which read [25] where r and r are the locations of the source and observation points, respectively. k = k 0 = ω √ 0 µ 0 is the wavenumber of free space herein. The term I is the identity dyad defined in global coordinates.
The region 2 only contains the internal fields generated by the IDTs placed outside of it, and the expressions of fields read where J 2 i1 , J 2 i2 and J 2 i3 are the three orthogonal electric point sources in the ith IDT associated to the region 2. N 2 is the number of IDTs placed in the region 1. Since region 2 is occupied by the uniaxial material, the two Green's functions, G 2 ee and G 2 me , in (6) are the uniaxial dyadic Green's functions, corresponding to the electric and magnetic fields radiated into the region 2 by an electric point source, which read [24], [25]

B. PLACEMENT OF MATCHING POINTS
The matching points should be distributed as uniform as possible in order to capture the behavior of fields on the physical surface of a scatterer with a smooth boundary. A simple and efficient way to place matching points is to make use of the Rao-Wilton-Glisson(RWG)-mesh [43] information which can be exported from commercial software package FEKO [44]. As shown in Fig. 4, the locations of nodes constituting a RWG-mesh are used to place the matching points. The parameter, TEL(triangle edge length) defined in FEKO, is used to control the density of a generated RWGmesh. In the MIDM, we firstly model the investigated object in FEKO and then run the mesh module to generate a RWGmesh. The matching points are followed to be placed with the help of nodes information of exported RWG-mesh, and a uniform placement can be finally achieved easily. The total number of matching points, N m , is determined by the TEL value defined in FEKO.

C. PLACEMENT OF IDTs
The IDTs are placed inside and outside of an object to construct the simulated fields outside and inside, respectively. A multiple-layered sources distribution scheme is proposed for the placement of IDTs in this paper. L virtual surfaces, which have the identical shape with the physical surface, are formed by scaling the physical surface to internal and external regions. L scale parameters are linearly spaced in the range 0.2 ∼ 0.995 for internal and 1.5 ∼ 2.5 for external regions based on empirical experiments. For the single-layered distribution, L = 1, the scale parameters are 0.2 and 2.0 for internal and external regions as found in [42]. A total number of IDTs (N s ) are unequally allocated to the L layers, and the number of IDTs in each layer is determined with respect to ratios between L scale parameters. The larger scale parameter corresponds to more IDTs. The placement of IDTs also makes use of the nodes information of a RWGmesh generated by FEKO. Once the number of IDTs for ith layer, N si , is obtained, we select N si points from N m nodes orderly, and transform those nodes to internal and external regions using the scale parameter of the ith layer. Finally, a multiple-layered distribution scheme of IDTs is realized. An illustrative example of three-layered IDTs distribution strategy is shown in Fig. 5.

D. BOUNDARY CONDITIONS
The connection between the fields in regions 1 and 2 is dictated by the boundary conditions of surface S shown in Fig. 1. Specifically, the tangential components of electric and magnetic fields must be continuous along the physical boundary S, which leads tô wheren is a unit vector normal to the closed smooth surface S as shown in Fig. 1. A matrix is then created by imposing the boundary condition at a number of matching points on S. Two tangential electric and magnetic fields are calculated at each matching point, therefore the formulated matrix has 4N m rows. On the other hand, 2N s IDTs in two regions are used to simulate the tangential fields, and therefore 6N s infinitesimal electric dipoles are deployed in the simulation resulting in a matrix with 6N s columns. The number of matching VOLUME 8, 2020 points (N m ) must therefore satisfy the inequality N m ≥ 1.5N s (9) in order to have a unique soloution for the unknown current coefficients.
Upon the application of a point-matching procedure, we will finally obtain a matrix expression of the type where X is a column vector containing the unknown dipole coefficients, and B is another column vector containing samples of incident tangential fields at the matching points.
[Q] is a matrix whose entries are obtained from the tangential fields of IDTs at matching points, and it could be rectangular or square depending on whether oversampling is used or not. If it is in a square form, a unique solution can be found, otherwise the smallest least-square error solution is pursued and known to be where Q is the transpose of [Q] and the asterisk denotes complex conjugate.

E. CONVERGENCE STUDY
In order to study the convergence of results, we make use of the error on the imposed tangential boundary conditions as a metric, whose definitions read where E bc and H bc are evaluated on a more dense points than the matching points selected on the physical surface S. The necessary numbers of sources and matching points in the MIDM are increased until E bc and H bc reach the desired level of accuracy.

III. SINGULARITIES IN USING DYADIC GREEN's FUNCTIONS
Two types of singularity issues usually appear when using the dyadic Green's functions. The first singularity issue is involved in both free space and uniaxial dyadic Green's functions when |R|, the distance between source and matching points, is approaching zero. A special treatment on this issue should be considered in the SIE-based solution whereas this type of singularity is completely avoided in our case because the sources are placed at a certain distance away from the matching points. The other singularity appears in the usage of the uniaxial dyadic Green's functions when the term |R c | vanishes, which means R is parallel to the optical axisĉ, and R c /|R c | becomes undefined as can be seen from (7). A special treatment on this issue has been proposed in [17] for the SIEbased methodology. However, that treatment is only valid for an electrically uniaxial object and should also valid for an magnetically uniaxial object according to the duality theorem. For a general uniaxial material, the proposed solution in [17] is not applicable. Moreover, up to authors' knowledge, the analytical solution for the second type of singularity of a general uniaxial material has not been reported so far. Instead of deriving complex analytical solutions, a quite simple approach to eliminate the second singularity issue by removing the problematic sources is proposed in our proposed method. Once the locations of matching points and sources are generated as specified in Sec. II-B and II-C, we implement a filtering strategy to find those locations of sources letting R to be parallel toĉ with a criterion |R c | 1e −5 . Those detected problematic sources are removed directly in the MIDM in order to avoid the second type of singularity issue. It was found that the number of problematic sources is much less than the total number of sources, therefore it is safe to remove them without scarifying the accuracy. This strategy has been found to be a quite efficient way to handle the second singularity issue in our numerical examples.

IV. NUMERICAL EXAMPLES AND OTHER DISCUSSIONS
Based on the numerical scheme described in the previous sections, a computer program has been implemented. The program computes the normalized bistatic radar cross section (RCS) (σ/λ 2 ) in xoz and yoz planes, defined as where E s is the total scattered electric field in the region 1, and λ is the incident wavelength. The MIDM procedure flowchart is shown in Fig. 6, and the initial value of TEL is defined regarding to the electric size of an investigated scatterer. For a 3D object, it is better to set the initial TEL value as 0.08 if the volume of the scatterer is within (2λ) 3 , and the initial TEL value is better to be set as 0.16 if the volume of the scatterer is between (2λ) 3 and (4λ) 3 based on the empirical experiences. When the volume is larger than (4λ) 3 , unfortunately, it is impossible for us give the initial value of TEL because the CST installed on the server we can access encounters a stagnation in the simulation and hard to generate simulation results. The step of the refine mesh operation in Fig. 6 is set as 0.01 to regenerate a RWG-mesh in FEKO. We use L = 4 for the multilayered distribution of sources in four illustrative examples, and the layer number factor will be analyzed and discussed later. The computed results are compared with the reference results obtained from commercial software package CST where 4 cells per wavelength for the model and the background is set in the mesh properties of CST.

A. FOUR ILLUSTRATIVE EXAMPLES
The first example is a both electrically and magnetically uniaxial sphere with radius r a = 0.5λ illuminated by a plane wave with an unit magnitude ofx polarized electric field and propagating along the z axis. The uniaxial medium parameters are 1 = 2 0 , 2 = 4 0 , µ 1 = 3µ 0 , µ 2 = 5µ 0 and the optical axisĉ is parallel to the z axis. The normalized bistatic RCS responses in xoz and yoz planes are shown in Fig. 7(a) and 7(b). The computed results are compared with simulation results obtained from commercial software CST [45] where the finite-element method (FEM) is applied. Good agreement can be observed between the MIDM using single-layered (L = 1) or four-layered (L = 4) distribution scheme and CST. To further quantify the performance of the proposed MIDM, the relative difference, as defined in (14), between the normalized bistatic RCS responses obtained from the MIDM with L = 1 or L = 4 sources distributions and CST is evaluated.

Error(dB) = |RCS(MIDM ) − RCS(CST )| (14)
The relative difference responses of the first numerical example are shown in Fig. 8(a) and 8(ab). Apart from a big difference appeared around θ = 90 • in yoz plane, the relative difference responses are small for the MIDM with L = 1 or L = 4 as observed in Fig. 8, which corresponds to the good agreement presented in Fig. 7. It is reasonable to have a big difference appeared around θ = 90 • in Fig. 8(b) because the relative difference is evaluated at the peak values of simulated results obtained from MIDM and CST. The second example is related to an electrically uniaxial capsule with a height h = λ and two end-capped hemispheres with a radius r a = 0.5λ, as shown in Fig. 9(a), illuminated by a plane wave with an unit magnitude ofx polarized electric field and propagating along −z axis. The uniaxial medium parameters are 1 = 5 0 , 2 = 9 0 , µ 1 = µ 2 = µ 0 and the orientation ofĉ is defined as θ c = 45 • and ϕ c = 90 • . The computed normalized RCS results in xoz and yoz planes of two IDTs distribution strategies are shown in Fig. 9, and are compared with those obtained from CST. For the singlelayer scheme, which is referred from [42], a disagreement on the normalized RCS computation is obvious to see, yet the proposed multilayered distribution strategy with L = 4 works well and excellent agreement is observed. This phenomenon is also indicated by the relative difference responses in xoz and yoz planes as shown in Fig. 10, where a small difference is observed for the L = 4 sources distribution scheme whereas VOLUME 8, 2020 a large difference is obtained for the L = 1 counterpart. Since the axisĉ is oriented in the yoz plane, the scattered field pattern in the yoz plane will be asymmetric as observed in Fig. 9(c).
Another different shape of scatterer is considered in the third example where six offset spheres merged. Six spheres are identical with a radius r a = 0.5λ, and are placed by offsetting one's center to six axes with a distance 0.5r a , as shown in Fig. 11(a). The uniaxial medium parameters are 1 = 2 0 , 2 = 8 0 , µ 1 = 3µ 0 , µ 2 = 9µ 0 , therefore a relative large ratio of anisotropy is considered, and the optical axisĉ is parallel to the z axis. This scatterer is illuminated by a plane wave with an unit magnitude ofx polarized electric field and propagating along z axis, as shown in Fig. 11(a). The normalized RCS results are computed in both xoz and yoz planes with single-and multiple-layered IDTs distribution schemes, and are compared with results obtained from CST as shown in Fig. 11. Once again, the single-layered scheme fails to have a good agreement on the normalized RCS simulation with that from commercial software package, whereas the proposed multilayered IDTs distribution scheme can. This also coincides with the results presented in Fig. 12 where a large difference is observed for the MIDM with L = 1 and a small difference is achieved when the multilayered sources distribution scheme is deployed in the MIDM.
The last example is about scattering from a uniaxial layer coated PEC sphere as shown in Fig. 13. In Comparison to previous three examples, three regions instead of two are presented in this scenario. The outermost region 1 is free space, the region 2 is a layer occupied with an anisotropic T i O 2 material, and the innermost region 3 is a PEC sphere with a radius r a . The medium parameters of T i O 2 are 1 = 5.913 0 , 2 = 7.197 0 , µ 1 = µ 2 = µ 0 as referred from [46], and r b = 2r a = 0.6λ. The incident plane wave has a unit magnitude ofx polarized electric field and propagating along z axis. The formulation is similar to previous examples. The only difference is that additional sources located within the region 3 are required to simulate the fields inside of the anisotropic layer, as shown in Fig. 13(b). It is noteworthy that the locations of sources in Fig. 13(b) are randomly depicted to illustrate the concept of MIDM, and the practical implementation of sources placements is referred to Sec. II-C. The simulated normalized RCS responses are presented in Fig. 14. Since the electric size of this example is relatively small, both L = 4 and L = 1 sources distributions in MIDM can generate normalized RCS results which have good agreement with that obtained from the CST. The relative difference responses of the two sources distributions are also small as presented in Fig. 15.

B. SINGULARITY ISSUE AND SIMULATION PERFORMANCES
The TEL value is set as 0.08 in FEKO to generate a RWGmesh for placements of matching points and IDTs in three numerical examples. The resulting numbers of matching points (N m ) and IDTs (N s ) for each example are displayed in Table 1. As discussed in Sec. III, two types of singularity issues are encountered in using the dyadic Green's functions as discussed in [17]. The first type of singularity issue is not existing in the proposed MIDM, but the second  type of singularity occurs. With the help of the freedom of sources placement in MIDM, a filtering strategy can be easily implemented to find those problematic sources as introduced in Sec. III. The number of problematic sources (N sp ) detected under the criterion |R c | 1e −5 in each numerical example is presented in Table 1. Specifically, in the fourth numerical example, the RWG-mesh generated by FEKO results in 854 and 218 matching points on physical boundaries S 1 and S 2 , respectively, and the numbers of IDTs are 569 for regions 1 and 2, and 145 for the innermost region 3. The numbers of problematic sources are 6 in region 1 and 4 in VOLUME 8, 2020  region 3 for single-layered sources distribution scheme, and are 4 in region 1 and 3 in region 3 for the fourth-layered one. Since the number of problematic sources is much less than the total number of sources, it is therefore safe to directly remove them in the MIDM resulting a singularity free numerical solution in simulating uniaxial materials. The simulation performance including the CPU time and the required memory of the proposed MIDM and commercial software CST is displayed in Table 2. The cost on the required memory is drastically reduced by using the MIDM in comparison to the FEM-based solver of CST. This is because a surface discretization method is deployed in MIDM whereas the volume discretization strategy is utilized in the FEM-based solver of CST. It is noteworthy to mention that the surface discretization-based method for simulating anisotropic materials has yet been incorporated in any commercial software so far. Except for the second example, the CPU time by using the proposed MIDM is also smaller than that used in CST, and the simulation performance has more distinct advantages by using the proposed MIDM if a relatively larger scatterer is encountered, as will be introduced in Sec. IV-C. It is noteworthy that our programs are written in MATLAB [47] and we authors are not professional programmer. It is expected that the performance of the MIDM can be improved if a compiled programming language is utilized.

C. LAYER NUMBERS AND BOUNDARY CONDITION ERROR RESPONSES
The four-layered sources distribution scheme has been proven to be feasible for the simulation of uniaxial objects with a smooth boundary and better than the single-layered one in previous examples. Yet L = 4 is a special case, it is therefore necessary to investigate the MIDM with different layer numbers in order to validate the feasibility of multiplelayered scheme. We then study a T i O 2 sphere with a radius r a = 2λ illuminated by a plane wave with an unit magnitude ofx polarized electric field and propagating along z axis. The TEL is set as 0.16 in FEKO to generate a  RWG-mesh for placements of matching points and sources, and resulting 2346 mathcing points and 1564 IDTs. The normalized RCS results are computed in both xoz and yoz planes under different layer numbers of sources, and are presented in Fig. 16. The number of problematic IDTs is much less the total number of IDTs for each different layer number scenario, and is omitted herein. The results generated by using L = 4, L = 8 and L = 30 are almost the same and have excellent agreement with that obtained from CST, whereas the signle-layered scheme fails to provide an accurate solution.
To further study this phenomenon, we turn to evaluate the boundary condition error response, which is a straightforward way to examine the accuracy of the solution. The tangential electric and magnetic field boundary condition errors, as defined in (12), in the xoz plane for single-layered and four-layered distribution schemes of IDTs are presented in Fig. 17. It can be observed that the multilayered scheme, L = 4, possesses a better response in E bc or H bc , which means a solution with a higher accuracy can be obtained if the multiple-layered strategy is deployed in the MIDM. By considering the inaccurate responses generated by the single-layered scheme as shown in Fig. 16, a boundary condition error criterion is therefore required to be set. The criterion ( E bc & H bc ) ≤ 0.1% is used in judging the acceptance of errors in MIDM procedure as shown in Fig. 6 to make sure the computed results are accurate. The simulations of the T i O 2 sphere problem were run on a large memory (384 GB) installed server with an Intel(R) Xeon(R) E5-2680@2.70 GHz, and the CPU time and required memory for the proposed MIDM are 8811 s and 1.31 GB, respectively, whereas for commercial software package CST are 82980 s and 170.20 GB, respectively. Obvioulsy, the simulation performance has been drastically improved by using the proposed MIDM.

V. CONCLUSIONS AND DISCUSSIONS
The MIDM, an efficient momentum solution, has been introduced, formulated and further utilized in the threedimensional scattering computation of uniaxial scatterers with a smooth boundary in this work. The uniaxial dyadic Green's functions have been deployed in the GMT-like method, namely the proposed MIDM, for the first time. A simple and efficient strategy to avoid the second type of singularity issue, which cannot be analytically solved so far for a general uniaxial medium, is proposed. The placements of matching points and sources, which play a key role in the MIDM, have been discussed and specified in detail by making use of the RWG-mesh generated by FEKO. The proposed multiple-layered distribution scheme of sources has been proven to be feasible and accurate in the scattering simulation from a relatively larger objects in comparison to the traditional single-layered counterpart. Several numerical examples are investigated under different scenarios, such as shapes of scatterers, electrical sizes of scatterers, and material characteristics, and the computed results of the proposed MIDM with a multiple-layered sources distribution scheme for each example have excellent agreement with simulated results obtained from commercial software package CST.
The cost on required memory has been drastically reduced by using the proposed MIDM, and the CPU time of simulation in the MIDM is also better than that in CST when a general uniaxial material or a relatively larger scatterer is considered. Only uniaxial materials have been investigated in this paper, yet the MIDM can also handle certain types of anisotropic materials, such as chiral, isotropic warm plasma and bi-isotropic materials, since their closed-form Green's functions are available [25]. Although the anisotropic materials with an arbitrary full tensor permittivity and permeability or with a general inhomogeneous configuration which can be handled by the volume discretization-based methods are highly impossible to be tackled by the surface discretizationbased methods, for example the proposed MIDM, the proposed technique has shown a significant advantage on the simulation performance in this work when a homogeneous general uniaxial material is considered, and this advantage is definitely inherited in electromagnetic computations of other anisotropic materials whose closed-form dyadic Green's functions are reachable by using the proposed MIDM.
The limitation of the proposed MIDM appears in simulating an object with sharp corners. The fast variation of fields near the corner requires more matching points and sources to be placed around the corner in order to better approximate the behavior of fields there, as discussed in [35], [48], [49] where a two-dimensional scenario was discussed. In this case, the proposed methodology may fail to apply. To overcome this drawback, an adaptive mesh strategy by applying a fine mesh around the corner and a standard mesh on other areas should be deployed, and a strategy to place more sources around the sharp corner is also required in the MDIM. Another alternative solution could be found in [36] where a RWG-mesh based testing method and a randomly distribution scheme of infinitesimal dipoles are deployed, and scatterers involving sharp corners has been successfully simulated by using the RAS method. The electromagnetic simulation of an anisotropic scatterer with sharp corners by using the MIDM would be our future work and will be presented in the future publications.